{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from osgeo import gdal, ogr, osr\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 栅格图像分区统计(shapefile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "shp = ogr.Open(\"region.shp\")\n",
    "lyr = shp.GetLayer()\n",
    "featList = range(lyr.GetFeatureCount())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "POLYGON ((12489962.8787558 4556776.48105855,12491479.4573327 4556401.54172477,12490959.9962259 4554623.3513856,12489296.6922533 4554793.23786456,12489962.8787558 4556776.48105855))\n",
      "POLYGON ((12501366.7011455 4557749.00490835,12503796.6238343 4556934.10816401,12503119.6185079 4555406.16833917,12501108.8587108 4555319.0750846,12501366.7011455 4557749.00490835))\n",
      "POLYGON ((12500818.5319626 4543599.76428432,12502980.4557383 4543538.13193089,12502803.0390906 4542222.9650258,12501158.2726934 4541486.72128036,12500394.8256903 4542027.76587492,12500818.5319626 4543599.76428432))\n",
      "POLYGON ((12492511.4886615 4547078.48045433,12494581.1903936 4546613.002575,12494304.3116179 4545295.43989119,12492557.3252735 4544656.61943763,12491938.1918222 4545502.77484844,12492511.4886615 4547078.48045433))\n",
      "POLYGON ((12503831.3893832 4549358.34805163,12505449.6495598 4548831.96994797,12504920.0680193 4547509.88349369,12503176.5317089 4546620.37619032,12501596.9154015 4547803.26025338,12502023.7418015 4549224.66121237,12503831.3893832 4549358.34805163))\n",
      "POLYGON ((12506011.9281801 4539757.93686265,12507472.6545453 4539531.56411803,12507501.7333681 4537867.54832073,12505161.5565255 4536716.52049997,12504613.8965942 4539279.00700054,12506011.9281801 4539757.93686265))\n",
      "POLYGON ((12487701.1406434 4538508.5118863,12489721.0457234 4537942.51907771,12489555.0446427 4536123.78556627,12488112.017598 4535440.24642181,12486743.0971965 4536119.43864344,12486421.3008937 4537172.13358731,12487701.1406434 4538508.5118863))\n"
     ]
    }
   ],
   "source": [
    "for FID in featList:\n",
    "    feat = lyr.GetFeature(FID)\n",
    "    geom = feat.GetGeometryRef()\n",
    "    print(geom.ExportToWkt())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'PROJCS[\"Mercator_1SP\",GEOGCS[\"GCS_unnamed_ellipse\",DATUM[\"unknown\",SPHEROID[\"Unknown\",6378137,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Mercator_1SP\"],PARAMETER[\"central_meridian\",0],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]]'"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raster = gdal.Open(\"dem1.tif\")\n",
    "raster.GetProjectionRef()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature = lyr.GetFeature(0)\n",
    "# feat.GetGeometryRef().ExportToWkt()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_origin: 12484418.383615404, y_origin: 4563153.762195187\n",
      "pixel_width: 42.975698519009455, pixel_height: -42.97569851900933\n"
     ]
    }
   ],
   "source": [
    "transform = raster.GetGeoTransform()\n",
    "x_origin = transform[0]\n",
    "y_origin = transform[3]\n",
    "pixel_width = transform[1]\n",
    "pixel_height = transform[5]\n",
    "\n",
    "print(f\"x_origin: {x_origin}, y_origin: {y_origin}\")\n",
    "print(f\"pixel_width: {pixel_width}, pixel_height: {pixel_height}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "source_sr = lyr.GetSpatialRef()\n",
    "target_sr = osr.SpatialReference()\n",
    "target_sr.ImportFromWkt(raster.GetProjectionRef())\n",
    "coord_trans = osr.CoordinateTransformation(source_sr, target_sr)\n",
    "geom = feature.GetGeometryRef()\n",
    "geom.Transform(coord_trans)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'POLYGON'"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "geom.GetGeometryName()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "ring = geom.GetGeometryRef(0)\n",
    "numpoints = ring.GetPointCount()\n",
    "points_x = []; points_y = []\n",
    "\n",
    "for p in range(numpoints):\n",
    "        lon, lat, z = ring.GetPoint(p)\n",
    "#         print(f\"point: {lon}, {lat}\")\n",
    "        points_x.append(lon)\n",
    "        points_y.append(lat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_min = min(points_x)\n",
    "x_max = max(points_x)\n",
    "y_min = min(points_y)\n",
    "y_max = max(points_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_count: 51, y_count: 51\n"
     ]
    }
   ],
   "source": [
    "x_offset = int((x_min - x_origin) / pixel_width)\n",
    "y_offset = int((y_origin - y_max) / pixel_width)\n",
    "x_count = int((x_max - x_min) / pixel_width) + 1\n",
    "y_count = int((y_max - y_min) / pixel_width) + 1\n",
    "\n",
    "print(f\"x_count: {x_count}, y_count: {y_count}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 在内存中创建一个栅格对象\n",
    "target_ds = gdal.GetDriverByName('MEM').Create('', x_count, y_count, 1, gdal.GDT_Byte)\n",
    "target_ds.SetGeoTransform((x_min, pixel_width, 0, y_max, 0, pixel_height))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 使对象栅格的投影坐标与原始栅格保持一致\n",
    "raster_srs = osr.SpatialReference()\n",
    "raster_srs.ImportFromWkt(raster.GetProjectionRef())\n",
    "target_ds.SetProjection(raster_srs.ExportToWkt())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 将区域多边形栅格化\n",
    "gdal.RasterizeLayer(target_ds, [1], lyr, burn_values=[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_offset: 113, y_offset: 148\n",
      "x_count: 51, y_count: 51\n",
      "(51, 51)\n",
      "[[1270. 1278. 1295. ... 1337. 1348. 1358.]\n",
      " [1269. 1276. 1293. ... 1337. 1347. 1359.]\n",
      " [1267. 1275. 1291. ... 1336. 1348. 1360.]\n",
      " ...\n",
      " [1216. 1221. 1236. ... 1332. 1335. 1340.]\n",
      " [1222. 1225. 1238. ... 1323. 1323. 1328.]\n",
      " [1226. 1229. 1237. ... 1313. 1313. 1316.]]\n"
     ]
    }
   ],
   "source": [
    "# 将栅格作为数组读入\n",
    "banddata_raster = raster.GetRasterBand(1)\n",
    "data_raster = banddata_raster.ReadAsArray(x_offset, y_offset, x_count, y_count).astype(np.float)\n",
    "\n",
    "print(f\"x_offset: {x_offset}, y_offset: {y_offset}\")\n",
    "print(f\"x_count: {x_count}, y_count: {y_count}\")\n",
    "\n",
    "print(f\"{np.shape(data_raster)}\")\n",
    "print(data_raster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fcfbb622630>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnV2MJNd13/+nqvpjPvaDu0vSNJcO5YAJ5IdYAhaKAuWBoayAUQzJD3JgxQgYhABfHECOHUhkAgRxEAPSi6WXwAEBCeaDY0qOHZAgbDgEIyIwEFBaSZRNmpZJC5S1Is2lRO7X7Ex3V9XJQ/fu9P3fM31rZmd7ZnXPD1js3Oqqe0993K4+554PUVU4jpMfxUEL4DjOweCT33EyxSe/42SKT37HyRSf/I6TKT75HSdTfPI7Tqb45HecTLmhyS8iD4rId0TkNRF5dL+Echzn5iN79fATkRLAXwH4CIBzAL4O4JOq+hc7HVOurWnv+InrbbW+enotjRNJHLaVdwCUtzXUDocwSQ5TGNeN9pEy3KegY8oiFqSUcFtF7X5RL9wfiM9fSbDG+M5PPQZ8+pU00T4sa8Ft0PUwbgT3UdIxVfRAACz6mDZsaS9o11pGfdT0MLZ0jSQaxd62qA9r9za6siF8nQd0/wFgIJPrf//tuRoX3uEHfoe+u+y0Ax8A8JqqfhcARORJAB8HsOPk7x0/gXt+5d9db9erxtW4fRQK2A9PvizDh6Ou4we5HtNpXQpvvozDa1PU8bXiOdX2QlmbFeMbZEATdX0ctFeHYfvE2tWoi2P9zaB9chDuc3r4brh/FfcxasPzbejBvtisRMfU7eIfgQV9G57qXYn2OVVdDtpHyvBchnMPKQCsFeG9BoCTRXg+R4rwmFNlPHFb+ub6fhOey1+O7wzab02OR338sF4P2lt0DXvGl531xbuoj9Z4UW02/aDNX5i39cLr8feGfxv18Xd756///W8+9oOFMgVjdd4z5m4A359rn5ttCxCRR0TkrIicbTY2bmA4x3H2kxuZ/NZPi+hVrqqPq+oZVT1Trq3dwHCO4+wnN/Kz/xyAe+bapwG8kTxq7utBDNWkbcNtwso30e/HP8eYyWr4HacVtcfGz36WjZutoXvStnoc/kQdleHlvrQ1NOWdp1+G58c64sTQXwf0U5n3OSG7/wXGP0dvp5/4APCTVaiSHC22Fva5SnICsY7Pb6d3mvh+v90OgvYb9W3h5/XRoM0/8QHg7fGRoD1q0lODVSGGf+ZPDCNX3Yb3pirC8+Of/db9nr83nZT9a8ftYl/m6wDuE5H3iEgfwC8BePoG+nMcZ4ns+c2vqrWI/FsAfwKgBPAlVX153yRzHOemciM/+6GqfwTgj/ZJFsdxloh7+DlOptzQm38vFHOGNF6fBQClNVo2AFYVOYEYjjLDYbhNaQ27YYNfERtRlA1+LGppGHssx5/5cenc6ib+7mUj0bgJZRtI6ORhrj+zsLyPcYxlSJpntQh9FG6vLkX73FGGa/+rhkNKCpb9soaPqLVW/ur4J4L2G7SOf35MBr9xbPC7OA6Nr+wIxA5KQNpBh30nLNnZIevkMDTGso+Gdb/7mDf4dXfa8ze/42SKT37HyRSf/I6TKcvV+RWhk4/lGk26NuvrqhTo0Iv1StajyrXQ2WRck/NNFfpgA0BLMQNsAxBDvy+qcBs7KFVVKHtVxhegT04e673Y/z3F1Tb0F2cf8y5wzADrmieL2FGI/fCPFGS/IBvPxFBPL5KsGxSU8/3JyeiYFzd+KmizTs/6/MYkdAqabguvWRfNmXX4JnpW0w5ra/3QlsIBRhw/UBrBUOO5dzgHcS3C3/yOkyk++R0nU3zyO06mLFfnlzAu3kzmwQkwSC9mPZkDXwBgUIV2gFEdnuYqqXxXq7iP0SQ8htfoLTjXACfvWOmHOvHRYRz4cmwQxsAf74Vt1gGvNrH+OqK18asUM74n6EnhBBnTbaH+WrThdZ2QPnqhjeW60KwG7bebcI3+3PgEmLdGYVDOFdLp3x2FfV4ZGdeM7EAtr9EbgVwM3++CfFDWBuH1AeLgoLUq3Ge9XBwcBQBX5wKbLF+CnfA3v+Nkik9+x8kUn/yOkylL1fm1AJqVbR2n7Rtr5b3FOft4XX+9H6+Dsx2A/aXZ5/piFee0uzqhpI/NYt93AKho3AG1jw5C/W29imW/a3gx7JP0xigvnLGu++4k1HE5T9xGnbYBDMrwOq+VoV59vIxzB47JH75PvgENybrRxrr35Ta8F+9Q4g32YQCAcRs+xhfHYR+XtsJxNkdGH5T3MYrtMOB1/N6Ans3V0MbDzwMQP5un+nFuxHms5KsX2tWFn++Ev/kdJ1N88jtOpvjkd5xM8cnvOJmyXCefQtGszRmwenGQQq+/OHBntRcaUVar2HHiWD80rK2U4TFdMrOywauL8wQbGi3Z5jlOBTqA2OmDA2rY4LVRx0azS7SNnV6uGga/KBiKDI0cYPS9/qmoDzbgcVEOTkxhGae4oMikTTssbdWhEZSduiZkrJ1wURcA7YQzOqcNvELPKgduleT0s9qLn4cjhtF3EWwABYCtYvv8rWpEO+FvfsfJFJ/8jpMpPvkdJ1N88jtOpvjkd5xM8cnvOJnik99xMmW56/ylojyyveZeGkk01lfDNfojlACB10qtBJdrZbgPV63lAgzFIA4w4kSKXajMjKQ7M2njNdmWkl5eoCCdjQ5BOryuv0nr4FwIpAvsG/F6FSfSvFSHiTI5KGlI94HX/YE4WUkUyGT4W3AFYR63C0oJWzmRrFWQhQvMcBAa+6RYzyr7dXD7ShNe04nECWvnzz9VSCQ4rvOejuP8WOGT33EyxSe/42TKUnX+olAMV8ZzbaNoReQfTXoj+emXRiEE1nu4CCXriL0itj0cLWK/+xSjPRTHYFjHn5BefGFMNoCJ5ese3tZx1E7r/HxVt8qwjwEl9wCAzSY8/8j2QvESlVF0kmMZ+N5Zvuu8LSqQyW0rUUeqMKtVxYOePe6Vi3SwXEBsw+BYjboIj5n347/GvJ3EE3g6jpPEJ7/jZEpy8ovIl0TkvIi8NLfthIg8KyKvzv6/7eaK6TjOftPlzf87AB6kbY8CeE5V7wPw3KztOM4tRNLgp6r/V0Tupc0fB3D/7O8nADwP4DNdBpw34LFDAxBnyd3i5BZFaOBqDUeRzSo0CvaL2DEi6NNwzmFHoLVOSRfCcTeNTLPzXBgbWYPJaWdMjkAXR6HTx9YkvoWcvKLuYOBLVSRqqvDzSySHRdsLjU9sEGQDIBA7SrEBazdOLLuCKvIIV4vucAzvE90Hw1h5mQx8bOAtJXweODENEDpPLcPJ505VfRMAZv/fscd+HMc5IG66wU9EHhGRsyJytrkU53p3HOdg2Ovkf0tE7gKA2f/nd9pRVR9X1TOqeqY8urrTbo7jLJm9Ovk8DeAhAJ+d/f9U1wPn9SLWiQCgYR2PFCk+xqrS26Nt/WKx4xB/DsTJNzk4yIIdUjhw5+Ik1JMvjWO9mQNo2GHn8iYFekyMa5jQ37tUo0lhBQdxABHr63yvtoxEqrwP214sndZynrlhTCWfoOvI152v0eVxnGy1JlvKFtlF+PwbI+nrSrHdrxUstRNdlvp+D8D/A/D3ReSciDyM6aT/iIi8CuAjs7bjOLcQXaz9n9zhow/vsyyO4ywR9/BznExZbpVeFYxGvbl2vE9LwTGs4RWUMIETKADTAKJgH9LxuZouBxMBcSKGiuwCA8NOEAVpUOINTrJh6YC8br81poIUW+TnUMc6sNL6M69Hd6GgghR8rxqjTy78wTYP/pyDVoB00QnWgfeC7OGVJ9Z1pjbbXzgYyoLtAjwK268sm8e8n8q+6vyO4/x44pPfcTLFJ7/jZIpPfsfJlKUb/Caj7SF1Ynz3sGGF7RulLv4ciDKsFFVo8BPqo9+PA0zGg9AQs9YLvROtLDScNTaqjktZd66O4sCfzVFo4KvJANhcpVtmJardg4EP/bAjNs5pB0eaERmvSjqmoesj7MEFO9hrHitDMmdaTpoErTGibYnMPgCKUThuPQrvzYi6sJyv2GDNmamqcvF9AcJncTdZp/3N7ziZ4pPfcTLFJ7/jZMpyK/Y0gF7ZHtJynLC2hTsk2gCUdHpWk9peqEdt9mPHkvE4vDTsCMROQBYN6bwb41DHZ/0eQGATAYB2TE4gpGcmr5eBVlb1Geon8WRw0hULKj6EsSx2vpoes1hjt86WbS3sbMWOMpPSCIai1yA/Q8VW/J5kUdstCuziMap4XM7wm3Jiqw27wXwf7uTjOE4Sn/yOkyk++R0nU5ar87eCcmNb7zFiYyCkKCVUQFgqjpaLfQVqSihkrZ2yaBu0Js/6vCkb9Tuq04k1W07OQb4QxSbp/Nb1oW2s45tqISej5DYH5Ri6J18TrsgUBWkZwvMxbDdgHRnolndjHisYjP1HhHR8TugJAOxyUGyG945tD03PGJfPj3V8ClJLJlrdRWITf/M7Tqb45HecTPHJ7ziZ4pPfcTLFJ7/jZIpPfsfJFJ/8jpMpPvkdJ1OW6uQjLdDb2PZqEKN4LhfGMXJmBKhxBuzEUq9S8MSEnS8MNxFy0OBAn9ZwpmAHFYb7MCvn0DZ26im5WLAleuQ7RBVoLccoTs8bOf2EH1sOSuw807aLg38shx2R8Bi+ppzsAogdkDhHSORslPIcM7Ac0jiRSsHZOzjJiJUBuGAHrGLh5xMjUcuVOcefZhfVmPzN7ziZ4pPfcTLFJ7/jZMrSdf5ya7ttFb5lnVYa0okoaKc1zoADd/YEqYWcSLNtLSVwsY7LdoIoiAcAyB5Rkh5ZjBMJThHr75F5wspf2SwO7GlrSsbJiVQBtBR0IgkbiIVEgTxk8zCCcriyc0rrTQXHAHGSFNs+xc/i4iSgbDcBAOUT5s85KMuwG9SDuY53kbzV3/yOkyk++R0nU3zyO06mLDmZB1BubjfLsZHMge0AnHyTFFYdxDoO2wHYF6Dldf3CUII5icKYdN4O35vCa7isa05i2Qsah5ObFKR7GnlIjHV+ksvyDeBtrFvSmnVrXbNIFjqXDuvrkQ8CHdMY/gWc9JJhHb+1isVwYlRewzeTzYbtkuwmwraXDglBokQrZBOw7m07n3jE1/kdx0nhk99xMiU5+UXkHhH5qoi8IiIvi8inZttPiMizIvLq7P/bbr64juPsF110/hrAr6vqN0XkCIBviMizAP41gOdU9bMi8iiARwF8ZlFHokC1ta0HFWNjn4QOpH3SgcxCnWGTXcyjohUdElpGOrAFF7dMFaAwdE8hOwCvJUexDobs0Ro16edirC1HmzgJKLclXbSD7SbR1ehQZJVlt45pyAdBqR3pwca9LK+SP8FV9q+Ix2X7VKSv0/NQcFwGYvuL+TzPf27M2Hkbl1HHdEeSb35VfVNVvzn7+zKAVwDcDeDjAJ6Y7fYEgF/oPqzjOAfNrnR+EbkXwPsBvADgTlV9E5h+QQC4Y7+Fcxzn5tF58ovIOoA/APCrqnppF8c9IiJnReRsvbmxFxkdx7kJdJr8ItLDdOL/rqr+4WzzWyJy1+zzuwCct45V1cdV9YyqnqlW1vZDZsdx9oGkwU+m1qEvAnhFVX9r7qOnATwE4LOz/59K9tWGgTvFJO30UQ8TgQ+G3Smq0tvFwJeCHTgspw82Eu4+Z0TkxMNGpcjJxzqXqJJxwppnHBQZ+DjQhzNZWONG5Wj44w5ydPFZ4XvDBj522DGcq8rNVEBVPCzfi6hP45iIxDMSOaxZVa7mn8VdPHNdrP0fAvCvAPy5iLw42/YfMJ30XxGRhwH8DYBf7D6s4zgHTXLyq+qfYucoyQ/vrziO4ywL9/BznExZbmCPAuVozsmnjhWUpp/Q8dnpx0pGyfsYiSeStOygQRVoDb2R9a095IlM6vgcTJLIBbED1kFkJ+GkmKw3Wxc+5aASOezsQXjjEE6cGTk58X2wivRusc6PhW3AeH65magWbR3DOn4ku6Hzz9+J3Txz/uZ3nEzxye84meKT33EyZbkJPBUoJ4t1/rYixSjS36lt6fx8VnvRi1k97ZDQMVVgpAs8TrTun1hbnnYSNvkaWap2QRs5aUoUs2JmEVkslhZ8b40CHIl7Zem01UZinZ/vpXGfeE1+PtEsEAakXYdtBxxP1CH2qe2F7ej5TtxLIAxc6+QXMcPf/I6TKT75HSdTfPI7Tqb45HecTFlyxR4NnHxMY13CoJfKzLtTvwuxDFUcC0RGIsvwZgX7BJ9zliLDIBT1y3KQkdTO3LJ7C2dsu1psALSKEyUrKvMxe3DysQyt7KATVXpuF7etfqtNruxsCRM2azbeRdWi032knNhMg19/TlY3+DmOk8Inv+Nkik9+x8kUn/yOkyk++R0nU3zyO06m+OR3nExZbjIPANJsr0la663KASYU6BP5ARiJOpLJO/hja3dOxsnJLMZGAk9a505WY7EK3fKaNPsXdEgKWUYld+kadkj6yadfcHKPaAygSATUIEo+2qFqb5sOqEol3hDyjWCfDQDR/e2UfJO76IWy1pSsuhnGxzQDSqJC9yZZFQhhcNBufFz8ze84meKT33EyxSe/42TKcnX+Fii3tpWrZiX+7kn59nfxdY4U1qgCBQ9q6e+Lk2pY9grT/3v+80Syxum4dAzrqx2SefDptlS0xLJXpGjZtZ/1exixC1GRDm52SILKthbj/Kuri/3wo3tnFb7gE+xAvRLKPzlCn6+GfTYrRvKa/uJxW44XMXbX+T5c53ccJ4VPfsfJFJ/8jpMpPvkdJ1OWm8wDgMxZo9j5AggzkQKGQa+Lo0wi826EVXyGs/Wys4lVOWWXmXZNwxMb/CLjVdow1RSLjZVdqvRGzkacRMW6ZmzwS7xazKQaCScnM/MuZdqN70OXBCgsSNi0qkU3AzL4rYfj1EfCgXQYDyz9xRlQ2k2+8MZOmtrBxt/8jpMpPvkdJ1N88jtOpiy5Sq+iGG/rPfWKlQWS2lGpmPQwkZ4orK9zsFDcaVT1pYOzSayfLnbQsZ1NFo/TRec3eu2wS9hvS8cUrK5aXSbsManAp+k+i/Vz65pxNZ2boePXRlAOJ+RkJx7W8Yth/NCUvcU6/yRyWEs4RnkCT8dxUvjkd5xMSU5+ERmKyNdE5Nsi8rKI/MZs+3tE5AUReVVEviwi/ZsvruM4+0UXnX8E4AFVvSIiPQB/KiJ/DODXAHxeVZ8Ukf8O4GEAv53qLFBZLPWFzQCL81KgMNbwW1Imq0TATVQ9Fmlds0s11EjHZ/3d1Hm5j3Cn+SrHO2Hp0uEOhvDRJtJfy3SACQfHJItlWMFRiXV9KwCnHC/W+aPCGMb5Rzo+reHzmj4ANPS6awecEYTkMmxLw+Hih7NHNoGmid/X4825bB77qfPrlCvXZJn9UwAPAPifs+1PAPiF7sM6jnPQdNL5RaQUkRcBnAfwLIC/BnBBVa+9y84BuHuHYx8RkbMicnY82dgPmR3H2Qc6TX5VbVT1fQBOA/gAgPdau+1w7OOqekZVz/R7a9YujuMcALuy9qvqBQDPA/gggOMi1z2+TwN4Y39FcxznZpI0+InI7QAmqnpBRFYA/ByAzwH4KoBPAHgSwEMAnkr1pSJo+9sWvaZvZPJJVOVNOY4AcZaZlOHJrBZM27pklImcSxLZfyzDXEEGvciYxQa/Pfj8WNWBk5YiNqJ1MNbFlY7p+nRxcmrS51tthQdxFl028FlBOpHBb5Xb8bjtYPHFl4qcfMr4ovXK8CL0q7BdNyW144e1nXMEkg4Zka/Rxdp/F4AnRKTE9JfCV1T1GRH5CwBPish/BfAtAF/sPKrjOAdOcvKr6p8BeL+x/buY6v+O49yCuIef42TKcgN7ijBjr6V7cYWeSI9kh40OQRtcwSXlOLTjtvmP9xCUww4qXewGEZH9wqick6g2E10PWM5FtEMiq661LaWv233sPiiHbUdsJ6qH4edW5ZwJ6/i0MDVfFed6PwPakHiVlobOzzr+gGwAR/vhzRpzthsAxVzUVRFFYO2Mv/kdJ1N88jtOpvjkd5xMWarOryKBfmatN7POy0klSq70avVBsRK8Vt6FVGUgU19PjJsMQEGs83LQUTkKdToxquW2URBOKLxlr7CqBy3CGjel00f2iQ46P5+/ZeOoh+FDMF5frOM3hq1pwjo+VdKxdH6t2JGB7i+t8/eq+MIPq/BBOtbfos/Dh2rcGDr/nMNIuYt1fn/zO06m+OR3nEzxye84meKT33EyxSe/42SKT37HyRSf/I6TKT75HSdTlhvYI6G/iZVFNVWhJkqy0ala7u4DTPZSHThy6qHEG8WIPh/FA7PzDF+jckQnbDi9CDn5RI5TlXHd2TEokovaHSoHpRKiWI5C8TMR7tMM4/dVKnBnsk5BOyuxrG2Prjs59TR9Q1baJsPw3lQctNOLPcNWyIlntQqjstaoXZfG+c95upWdyhNN8Te/42SKT37HyRSf/I6TKcvV+Qkr4CZVbSau4BLvwzpuqmqrFSzCcFJIS+ePxuVknJsNfW5F2HCwCOnvm6GOqIXx/c1BSaQnFobeyOMwXKGIq+nupQ+zchDFrbSUqCOqHIQ4qQZX1+HKOlFADgClmRDbAIxjehS40w91+gG11/pxlpWjFMjDOv7RajNoj4wIo81qe1uRLNe0jb/5HSdTfPI7Tqb45HecTFmqzi+totrc1pMiPRqAXI2PCUgliQRQTLgqB/XJxSM66Pxd9NUo+WRDiTdojb4YGRVaee2bx5mEeqSYOj8VrShJkba+8mkc05aQ6EN7ZjWQbcjWoIbhhMflxCRNP77unNAlSvASVemNRUslb1Gjwi5fg4qSd6wOQv2d1/QBYKUMt7GOv16FCTyrNjZyjeYysRS7qOLib37HyRSf/I6TKT75HSdTlqzzA9VVY217fp+aE1TS57Q2zno1gLiwRcP+8CxX2h860oGtr83kuGQTGCWqawAAj8t9WkQ6fr34cyDSxwXhOHz+Okjo9wDaio7htuEX0AxI56c1e8tOxLEKVuHVmwIn7KR2SQ9vP6pAA1TkdMI6/noZ+gEMuPorgHruurpvv+M4SXzyO06m+OR3nEzxye84mbJcg1/dovfO1eQ+AWzQI4OX1IYBjB1lUk48HQx+UpGBywpKYTHomC7ORBFs4LPON4FQtIya0VAcuUNtcrZp+/Gj06yE29g4xwY+NZKKsINOzQZA64ll0RPVki1DI1uWo0pBRkCRFWS0iMK47gMSjp10jpWh08+kWGxoZQPiQnk67+k4zo8VnSe/iJQi8i0ReWbWfo+IvCAir4rIl0Wkn+rDcZzDw27e/J8C8Mpc+3MAPq+q9wF4F8DD+ymY4zg3l046v4icBvDPAfwmgF8TEQHwAIB/OdvlCQD/GcBvL+yoblD86NJ229K1E04tWoc6klrVYhN6EbS7XnQNaehSWTr/YPGPHx1yVglDdu6X9+HAHsMGoD2SlWwPVgCO9llfp2PIQWeyHj867JDTJHT+KCkqgDYRlMN97gUr6Wuc4IWPMXR+K0JojqpoF7aB2A6wWoZOPn1y6rmjugRmrdg+pmc4Ae1E1zf/FwB8Gts+bCcBXFDVayOdA3B351EdxzlwkpNfRH4ewHlV/cb8ZmNX05QtIo+IyFkROTtuN61dHMc5ALr87P8QgI+JyEcBDAEcxfSXwHERqWZv/9MA3rAOVtXHATwOAMf6d+5hrctxnJtBcvKr6mMAHgMAEbkfwL9X1V8Wkd8H8AkATwJ4CMBTyb7qCZq3zi/cRyoSqUcJCydh8gPpG3o25zhkG4Bwpobd2wCwMow26SAcuDkS7sO6NienBOKiFewbUF0c0eex7DxOM6yoHev8qUSZHDxTD+Mff5NVDsIJP+e4FisAJ9LH2QRiHMPjsC9AQjWfjbv4unfJkVGQTi/kO2Al1zxahYE7q0UY7PWT1bvh/kW4PwD8RHnx+t/DwkgQswM3ss7/GUyNf69hagP44g305TjOktmVh5+qPg/g+dnf3wXwgf0XyXGcZeAefo6TKT75HSdTlluxR0Mnnci4B0QGPuHqMsMj4ecccANEzjYcYJN0pEHsPMN96CqViQHQrIXb6rWK2mTwMwJb2MDFVX+i3Y3sxc0KjdNjY51V6TYRhEOXkKvkAEDNBj9OKMTGPMvHKRG3ZBn8WJbI4Mf2XaNiT+oYO+Px7jL3rJVx5iY28N1OTjx3lFeC9qkyNujNX+YBp5NagL/5HSdTfPI7Tqb45HecTPHJ7ziZ4pPfcTLFJ7/jZIpPfsfJlOUm8CxLlEePbbeHRnDM+mq4oR+u+0dBKyscxWNUiulRm9b5y1G8uBxV2KVKQe0wHpfX9cdHQlnrlcXJLizKMQWH1FySNj6mpUq2nAST5QCAmm5FlHwzkWTD6oPXyqOYHWudP7FMbQXpNCthRxzok6rAa45DvgCWbwBoW1WGwg+qMLHGGlXjAYBjVZjQ9mRiXf/2MnawGMj2Cffk3ejznfA3v+Nkik9+x8kUn/yOkynL9e0vC8ixo9ebOoz1l/rEWtBuB6w3c0KMdNXWyIee61GM4stQjkP9rdwM243hH8+yTFbDfRrKO9IaeUhSPvWTcfr7uulzUg2Sa824ZiQLyxrpyZbuPeTCFwvF7JbMI9oh3lSv0rg9apdc6jkxhoWh85fDUFjW8bniLhfoAMLkm12YaHyBirl3uHbJOnL9OMdxssQnv+Nkik9+x8kUn/yOkynLNfgVBXR12xOkXY0tXuPbwm31CjmoDBPGPMSJGRoj0+w85Sg2klRb4TEVGfgsZ5RmwAkxws+V5LJkjw1rZKybpL+vuQ825rETDBAb+KJkHdynlRCDjYQlO9/sg+HNQId0M6qwXVCbkzcDQGtU5JmnrIwKu8PQAWe9vzvjHQC8XYfJaYayOPvuySKucn2s2K6HMd5FJmp/8ztOpvjkd5xM8cnvOJmyVJ1fC0E7l+TSCsqpyTFmQjp/wwEoxhmwjs8BJ6xrVpuG3YBESzqfIHbISQaYGLKnqtTWFPdk+XRwtZ3IBmDo/HwdWT9XrorD+juAtr/Y2QY90keL3VdvE+OYivotq/Bm9fvpyrVNQzYduu6DXqyLD3thv6tVmIxzhYJyRsbDepWMK9/g4gVBAAAGYklEQVTTU0F7ow0/f9twCjoxFww00gvR5zvhb37HyRSf/I6TKT75HSdTlrzOL2gH20O2PSM4hta+o/VnClphGwAQr1FHgR/RmnWs85c0TknFUa1EFJxoIrIBRAFHRh90TB3p2jxI3AcXx2A7QpSYw9gW+SRw0Y5VYz2Z1v6lHwpSkG7OVW2BWNfmSrcWK4NQt+6Rzr9q6OspeNxBGdsNWMc/OYjX4AO5DMPRW+OjQbuifX4wOh602Y4AACeqjet/X21/sFCGefzN7ziZ4pPfcTLFJ7/jZMpy1/lFgmQcps7P+mgi+aK5Zp1ImhEngTQKdbZkWyA7QlR00oB1/EiP7rDezuv8fC6W/0FLqjTHIUTXGPG6fbTuz/7znCADaR2/Il28MNbsUzp+WcZ2gpSOz/o6J90AgLoNH6z1XrieXhnBHKx/c7KOgo65VK9EfYzoZlyZhA9aZdhFmB9W69f/3rIeqh3wN7/jZIpPfsfJlE4/+0XkdQCXMU29XqvqGRE5AeDLAO4F8DqAf6Gq3ZOGO45zoOzmzf9PVPV9qnpm1n4UwHOqeh+A52Ztx3FuEW7E4PdxAPfP/n4CwPMAPrPwiAJo+9vfN5HDCowkGWT/iYJjjBwMbbk4e2vUR2sk1eCEIJyoo5PBj/pIZcQFoq/jVEZcK6mI1Iv3YWMmYATuDMKDooo1HKRjwKM0VG1IDeMdG/S6OPkoPQQTMt4NrWo7xNF+6MXFxjzLuWalCJ182Di32YQXtTUeVjbwbURtI8UzMW/AtIKHdqLrm18B/G8R+YaIPDLbdqeqvgkAs//v6Dyq4zgHTteviQ+p6hsicgeAZ0XkL7sOMPuyeAQABivHE3s7jrMsOr35VfWN2f/nAfwvAB8A8JaI3AUAs//P73Ds46p6RlXP9Ppr1i6O4xwAorpYHxKRNQCFql6e/f0sgP8C4MMAfqSqnxWRRwGcUNVPJ/p6G8D3AJwC8MP9OIElcKvIeqvICdw6st4qcgLbsv4dVb29ywFdJv9PY/q2B6Zqwv9Q1d8UkZMAvgLgpwD8DYBfVNV3Og0qcnZu1eBQc6vIeqvICdw6st4qcgJ7kzWp86vqdwH8rLH9R5i+/R3HuQVxDz/HyZSDmvyPH9C4e+FWkfVWkRO4dWS9VeQE9iBrUud3HOfHE//Z7ziZstTJLyIPish3ROS12fLgoUFEviQi50XkpbltJ0TkWRF5dfb/bQcp4zVE5B4R+aqIvCIiL4vIp2bbD5W8IjIUka+JyLdncv7GbPt7ROSFmZxfFpG0D+uSEJFSRL4lIs/M2odSVhF5XUT+XEReFJGzs227uv9Lm/wiUgL4bwD+GYCfAfBJEfmZZY3fgd8B8CBtO6zBSzWAX1fV9wL4IIBfmV3LwybvCMADqvqzAN4H4EER+SCAzwH4/EzOdwE8fIAyMp8C8Mpc+zDLemPBdqq6lH8A/hGAP5lrPwbgsWWN31HGewG8NNf+DoC7Zn/fBeA7By3jDnI/BeAjh1leAKsAvgngH2LqjFJZz8UBy3h6NmkeAPAMprFJh1XW1wGcom27uv/L/Nl/N4Dvz7XPzbYdZg598JKI3Avg/QBewCGUd/Yz+kVM3b+fBfDXAC6o6rVQtMP0HHwBwKcBXAvPO4nDK+sNB9stM4efVQDdlxpuABFZB/AHAH5VVS8JJ70/BKhqA+B9InIcU0/R91q7LVeqGBH5eQDnVfUbInL/tc3Grgcu64w9B9tdY5lv/nMA7plrnwbwxhLH3wudgpcOAhHpYTrxf1dV/3C2+dDKq6oXMM358EEAx0Xk2ovnsDwHHwLwsVnWqicx/en/BRxOWaE3EGx3jWVO/q8DuG9mPe0D+CUATy9x/L3wNICHZn8/hKlufeDI9BX/RQCvqOpvzX10qOQVkdtnb3yIyAqAn8PUmPZVAJ+Y7XbgcgKAqj6mqqdV9V5Mn83/o6q/jEMoq4isiciRa38D+KcAXsJu7/+SjRQfBfBXmOp9//GgjSYk2+8BeBPABNNfKQ9jqvM9B+DV2f8nDlrOmaz/GNOfn38G4MXZv48eNnkB/AMA35rJ+RKA/zTb/tMAvgbgNQC/D2Bw0NeU5L4fwDOHVdaZTN+e/Xv52lza7f13Dz/HyRT38HOcTPHJ7ziZ4pPfcTLFJ7/jZIpPfsfJFJ/8jpMpPvkdJ1N88jtOpvx/ESZTXDhRx2IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(data_raster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No data value: None\n",
      "[[0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " ...\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "band_mask = target_ds.GetRasterBand(1)\n",
    "data_mask = band_mask.ReadAsArray(0, 0, x_count, y_count).astype(np.float)\n",
    "\n",
    "print(f\"No data value: {band_mask.GetNoDataValue()}\")\n",
    "print(data_mask)\n",
    "# np.shape(data_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fcfbb57f320>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAADN5JREFUeJzt3V2MXdV5xvH/E+OPQBKBCSAX05pIpCIXBKQRoaIXFEJDaRS4IFVQVLmSJd+0ElFSgWmlSqlaidwEbqpWVkHxRRogHxUIpXItF1RFqgzmMxAXTBBNLCPcUlBIUV1M3l6cbXUyzHTOzOzzNev/k0bn7D17e7/2nGfWXuusdZyqQlJ7PjDpAiRNhuGXGmX4pUYZfqlRhl9qlOGXGmX4pUYZfqlRawp/khuSvJjk5SR7+ipK0uhltTP8kmwAXgKuB44BTwC3VtWPljpnUzbXFs5a1fXWg49f9s6Kz3npuTNHUInWq//mv/ifOplhjj1jDde5Eni5ql4BSHI/cBOwZPi3cBafynVruORs27//mRWf85lfuXwElWi9OlQHhz52Lbf9FwI/nbd9rNv3S5LsTnI4yeF3ObmGy0nq01rCv9itxfv6EFW1t6rmqmpuI5vXcDlJfVrLbf8x4KJ529uB42srRwvtP77yrsJCdh20mLW0/E8AlyS5OMkm4AvAw/2UJWnUVt3yV9WpJH8E7Ac2APdV1Qu9VSZppNZy209VfR/4fk+1SBojZ/hJjVpTy6/Z4KChFmPLLzXK8EuNMvxSo+zzj1Affe1p4bjB+mPLLzXK8EuNMvxSo+zza2wcN5gutvxSowy/1CjDLzXKPr9miuMG/bHllxpl+KVGGX6pUYZfapQDfj1aTwt51rO+fk6zPnBoyy81yvBLjTL8UqPs80urNOsTjmz5pUYZfqlRhl9qlH1+aYImOW5gyy81yvBLjTL8UqPs86+S8/g1Lea/Fq/8zDtDn2fLLzXK8EuNWjb8Se5LciLJ8/P2bU1yIMnR7vGc0ZYpqW/DtPzfAG5YsG8PcLCqLgEOdtuSZsiyA35V9c9JdizYfRNwTfd8H/AYcEePdUka0vxJPi/VG0Oft9o+/wVV9RpA93j+Kv8cSRMy8rf6kuwGdgNs4cxRX07SkFbb8r+eZBtA93hiqQOram9VzVXV3EY2r/Jykvq22vA/DOzsnu8EHuqnHEnjMsxbfd8C/gX49STHkuwC7gKuT3IUuL7bljRDhhntv3WJb13Xcy2SxsgZflKjXNgzJBfyaL2x5ZcaZfilRhl+qVGGX2qUA37SDOnzf/ix5ZcaZfilRhl+qVGGX2qU4ZcaZfilRhl+qVG+z78EF/JovbPllxpl+KVGGX6pUYZfapThlxpl+KVGGX6pUYZfapSTfKQp1ueHdyxkyy81yvBLjTL8UqPs8+MiHrXJll9qlOGXGmX4pUYZfqlRhl9qlOGXGrVs+JNclOTRJEeSvJDktm7/1iQHkhztHs8ZfbmS+jLM+/yngK9U1VNJPgw8meQA8AfAwaq6K8keYA9wx+hKlda/Uc7lX2jZlr+qXquqp7rnbwNHgAuBm4B93WH7gJtHVaSk/q2oz59kB3AFcAi4oKpeg8EvCOD8vouTNDpDhz/Jh4DvAl+qqp+t4LzdSQ4nOfwuJ1dTo6QRGCr8STYyCP43q+p73e7Xk2zrvr8NOLHYuVW1t6rmqmpuI5v7qFlSD5Yd8EsS4F7gSFV9fd63HgZ2And1jw+NpMIRcCGPNNxo/9XA7wM/THI6NX/CIPQPJtkF/AT4/GhKlDQKy4a/qn4AZIlvX9dvOZLGxRl+UqMMv9Qowy81yvBLjTL8UqP8AE9pgsa5kGchW36pUYZfapThlxpl+KVGNTHg50Ie6f1s+aVGGX6pUYZfapThlxpl+KVGGX6pUYZfapThlxpl+KVGGX6pUYZfatS6nNvvXH5No0l+cMdibPmlRhl+qVGGX2qU4ZcaZfilRhl+qVGGX2qU4ZcaZfilRhl+qVGGX2rUsuFPsiXJ40meTfJCkq92+y9OcijJ0SQPJNk0+nIl9WWYhT0ngWur6udJNgI/SPIPwJeBu6vq/iR/A+wC/nqEtS7KRTyaVtO2kGehZVv+Gvh5t7mx+yrgWuA73f59wM0jqVDSSAzV50+yIckzwAngAPBj4K2qOtUdcgy4cIlzdyc5nOTwu5zso2ZJPRgq/FX1XlVdDmwHrgQuXeywJc7dW1VzVTW3kc2rr1RSr1Y02l9VbwGPAVcBZyc5PWawHTjeb2mSRmmY0f7zkpzdPf8g8GngCPAocEt32E7goVEVKal/w4z2bwP2JdnA4JfFg1X1SJIfAfcn+QvgaeDeEdYpqWfLhr+qngOuWGT/Kwz6/5JmkDP8pEYZfqlRhl9qlOGXGmX4pUbN3P/Y40IeqR+2/FKjDL/UKMMvNWrm+vzStJr2D+9YyJZfapThlxpl+KVGGX6pUYZfapThlxpl+KVGGX6pUVM/yceFPNJo2PJLjTL8UqMMv9Soqe/zS9Nq1hbyLGTLLzXK8EuNMvxSowy/1CjDLzXK8EuNMvxSowy/1Kipm+TjQh5pPGz5pUYNHf4kG5I8neSRbvviJIeSHE3yQJJNoytTUt9W0vLfBhyZt/014O6qugR4E9jVZ2GSRmuo8CfZDvwu8LfddoBrge90h+wDbh5FgZJGY9iW/x7gduAX3fa5wFtVdarbPgZc2HNtkkZo2fAn+SxwoqqenL97kUNrifN3Jzmc5PC7nFxlmZL6NsxbfVcDn0tyI7AF+AiDO4Gzk5zRtf7bgeOLnVxVe4G9AB/J1kV/QUgav2XDX1V3AncCJLkG+OOq+mKSbwO3APcDO4GHRlinNFGz/sEdi1nL+/x3AF9O8jKDMYB7+ylJ0jisaIZfVT0GPNY9fwW4sv+SJI2DM/ykRhl+qVETXdjjIh5pcmz5pUYZfqlRhl9qlOGXGmX4pUYZfqlRhl9q1NR9gKc0DdbjQp6FbPmlRhl+qVGGX2rUWPv8H7/sHfbvdz6/NA1s+aVGGX6pUYZfapThlxo11gG/l547c8WTJ/zAD2k0bPmlRhl+qVGGX2rU1C/smaUFFo5PzK5Zep31xZZfapThlxpl+KVGTX2ff5bMSr/RsQmBLb/ULMMvNcrwS42yz9+gWRmbAMcnRsmWX2qU4ZcaNdRtf5JXgbeB94BTVTWXZCvwALADeBX4vap6czRlSurbSlr+36qqy6tqrtveAxysqkuAg922pBmxlgG/m4Bruuf7gMeAO9ZYj/RLZmlwctYM2/IX8I9Jnkyyu9t3QVW9BtA9nj+KAiWNxrAt/9VVdTzJ+cCBJP867AW6Xxa7AbZw5ipKlDQKQ7X8VXW8ezwB/D1wJfB6km0A3eOJJc7dW1VzVTW3kc39VC1pzVJV//8ByVnAB6rq7e75AeDPgeuAN6rqriR7gK1Vdfsyf9a/A/8GfBT4jz7+AmMwK7XOSp0wO7XOSp3wf7X+WlWdN8wJw4T/Ywxaexh0E/6uqv4yybnAg8CvAj8BPl9V/znURZPD8941mGqzUuus1AmzU+us1Amrq3XZPn9VvQJ8cpH9bzBo/SXNIGf4SY2aVPj3Tui6qzErtc5KnTA7tc5KnbCKWpft80tan7ztlxo11vAnuSHJi0le7t4enBpJ7ktyIsnz8/ZtTXIgydHu8ZxJ1nhakouSPJrkSJIXktzW7Z+qepNsSfJ4kme7Or/a7b84yaGuzgeSbJpknfMl2ZDk6SSPdNtTWWuSV5P8MMkzSQ53+1b08x9b+JNsAP4K+B3gE8CtST4xrusP4RvADQv2TevipVPAV6rqUuAq4A+7f8tpq/ckcG1VfRK4HLghyVXA14C7uzrfBHZNsMaFbgOOzNue5lrXttiuqsbyBfwGsH/e9p3AneO6/pA17gCen7f9IrCte74NeHHSNS5R90PA9dNcL3Am8BTwKQaTUc5Y7HUx4Rq3d6G5FngEyBTX+irw0QX7VvTzH+dt/4XAT+dtH+v2TbOpX7yUZAdwBXCIKay3u41+hsH07wPAj4G3qupUd8g0vQ7uAW4HftFtn8v01rrmxXbj/Ay/LLLPtxrWIMmHgO8CX6qqnyWL/RNPVlW9B1ye5GwGM0UvXeyw8Vb1fkk+C5yoqieTXHN69yKHTrzWzqoX2502zpb/GHDRvO3twPExXn81hlq8NAlJNjII/jer6nvd7qmtt6reYvCZD1cBZyc53fBMy+vgauBz3adW3c/g1v8eprNWag2L7U4bZ/ifAC7pRk83AV8AHh7j9VfjYWBn93wng771xGXQxN8LHKmqr8/71lTVm+S8rsUnyQeBTzMYTHsUuKU7bOJ1AlTVnVW1vap2MHht/lNVfZEprDXJWUk+fPo58NvA86z05z/mQYobgZcY9Pv+dNKDJgtq+xbwGvAug7uUXQz6fAeBo93j1knX2dX6mwxuP58Dnum+bpy2eoHLgKe7Op8H/qzb/zHgceBl4NvA5kn/my6o+xrgkWmttavp2e7rhdNZWunP3xl+UqOc4Sc1yvBLjTL8UqMMv9Qowy81yvBLjTL8UqMMv9So/wXWdEZ77+ZIoQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(data_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fcfbb556d68>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHd5JREFUeJztnXuM5WV5x7/Puc3Mzi4uC4MSdnEhQYsmXuoGtTSRojR49w80WlvXSEJTawJeIlDTpmo12lbBNI2GBiNNrXhtoMRLCZU0GoPuCii44iIF3ULYEXZZYG7n8vSPc1bn/b7vzO83Z86ZOWff7yeZzLy/63N+Z57znud5n4u5O4QQ+VHZbAGEEJuDlF+ITJHyC5EpUn4hMkXKL0SmSPmFyBQpvxCZIuUXIlPWpfxmdrGZ3Wdm95vZVYMSSggxfKzfCD8zqwL4BYCLABwC8CMAb3P3n610zqmnnuq7d+/u634nAgeO/V8wrlfa0TFV6wTjs6Z3D1MkcYKxf//+37j7TJlja+u4z3kA7nf3BwDAzG4E8EYAKyr/7t27sW/fvnXccrw579tXB+PTpp+KjtnemAvG//bS64cqkzixMLOHyh67nq/9ZwD49bLxod42FuYyM9tnZvtmZ2fXcTshxCBZj/JbYltkQ7j7de6+x933zMyU+jYihNgA1vO1/xCAXcvGOwE8vD5xTmzmlurB+PHKlsJzLv3RO4Px87eGfoPt1dBMAIAKQr/BO5/zg5ISipxYz8z/IwDnmNlZZtYA8FYANw9GLCHEsOl75nf3lpm9B8B3AFQBfN7d7x2YZEKIobKer/1w928C+OaAZBFCbCCK8BMiU9Y184u10W6Hn7WtTvzZy9sWO9VgPGnNYMzOPQBoWBg89IVfvLxQtk7BPPCu53y/8BpivNDML0SmSPmFyBQpvxCZIpt/iDz36x8JxpON0BavVRL2OiX7nFRbXPN9j3WmgnHTqyscuTLbKvPB+Kv3vyQY//5EHM81Uw3/nbZWJoNx5VkH1yyHGB6a+YXIFCm/EJki5RciU2TzD5GpiaVgvKURrtHvmIqTck6ZeDocN8Kc/4qFiZNznYnoGgseJhDNteNjCin4z1jweN5Y8NBfUemEr7/58M5g/EQnLmbyeDuU/XB7azC++OwVy0WINaKZX4hMkfILkSlSfiEyRTb/EDlta2ivb58I186fUQ/HALBr8kgwrlOcPq/Zp9bwj7Smg/FcuxGMj7XC9fcUE5WTgvFJtYVgfEotrj+46E+E17BWKCvNNXMUjwAAxzqhbI+Tzf/pAxdF59z7VFg97tGFbcH4sfmwaMoP/viT0TVyRDO/EJki5RciU6T8QmSKlF+ITJHDb4C85Qd/Hox3Ta+elLM94fCrUPVzdvg91Q4dYk+2Y+fd0Wbo4DrWCoN85lqhAxAAWlQ0pFENnXVP18JrPNiIy7DzMdOV8PW3KTCo6fG/35Pk8FvoFAcszVNg0GIrvO5CMxxzwhUANJfCY9pLxfPiQ3vHu0OdZn4hMkXKL0SmSPmFyBTZ/EL0ydlf+ngw3rY19OFsmwx9Hs+aPhZd42t/8LnBC1YSzfxCZIqUX4hMkfILkSmy+QfIrqkjxQcto+Nxl/MOdT5/nJJ0OCnn6Va87v0kbXuqGY4X2yXedjpnoR6upd8396zolCONUFaOUZiohMVMOomCIBULi5ou0jo/P5/UOWZRp/hCOk2SpUXjanzNTjuUpUIFWbfUw2Im2xLFWN9/11uC8an1MGGKnxk3bQGAd//e7dG2MmjmFyJTpPxCZIqUX4hMkc3fJx+55/XRtolKPXHk2jhCcfls8/5miQp1JOL0n26G2zjWfald3MSj3Qnt2QW6xmQ1tj3nqWjISbVw3ZublExW4msUFS/hWH8AWOqEsjn5Unjcasav3+n1Rv1PY1dD5FuoVsIxF1ttJnwc7PfhXA1+/YsWv/7P3feK3/595vO3viQ6YAU08wuRKVJ+ITKlUPnN7PNmdtjM7lm2bYeZ3WpmB3u/Tx6umEKIQVNm5v8CgItp21UAbnP3cwDc1hsLIcaIQoefu/+Pme2mzW8EcEHv7xsA3A7gygHKNfKwYwpA9FGaKrSxnKPNuHotB+0sUZGNxxZChx874oC4eEWzjIOvvfo80OmE+5+ox7JH55Azq1bhoJ/YeTVRCYuI8DX66ThcCgrYAQdgJd5up2fGTtI2PTN2TAJx4ZXFTugE5QCmqYSj9ernf/O3f//FvbY/ljRNvzb/M939EQDo/T6tz+sIITaJoTv8zOwyM9tnZvtmZ2eHfTshREn6Vf5Hzex0AOj9PrzSge5+nbvvcfc9MzNx3TchxObQb5DPzQD2AvhE7/dNA5NoTEjZnryN7dXHKUDniWbsE+CgnflWaBc/MR+es9iM38JWa3W7uFNg35eB5QKAKgXxxDZ/uL9B9j0ATJBfgG3eVDJQKkFq3XCQT+KRcmAQ+0X4GVWX4ve7Uw+v8TQFStXo9aeKvvZLmaW+LwH4AYDnmtkhM7sUXaW/yMwOArioNxZCjBFlvP1vW2HXKwcsixBiA1GEnxCZosSekrx7/58G46lqvM7NtuexVngMF9k4shiu8QLxuv3cYmgDPj0fjjuJNXwne5XXo8tQbYS2NyfHtFOFSGgbr2u3OmGiS6caXyO1Fr4ctoH7wSp9XIPjABKw/4WLeyTPocIq/PrYj9JK+Dz6RTO/EJki5RciU6T8QmSKlF+ITJHDryRcQadRiavm1slZww6+YxTk8fRSogoPOfgWFylJZ56Ca1KOqBLOKcYmQtm5Vi07/FKw84rP4co3nLQExIFADHcxBuKKvkVOMUu9lKjiL8meev1L5OBcDN8bfr0cBATETsF6NRzz8/j2Kz4Ty9EnmvmFyBQpvxCZIuUXIlNk85fk0bltwXhbI+6+wnCV3KeWQh8A2/cAsEBBPO0lsosXw7E1127fez22m6PqtQU2fpkCIXwMV7PlAJbUMUwqyCeyrUn2Ri0MWFqsxZU5lrgjD0+Li/HzYFHa81RpF+F72azG963SfdkHUK2uP6hpJTTzC5EpUn4hMkXKL0SmyOZfgbO++PFgvH17aK9zccYUbHvOU+LH0lL8+NvUHdape2xlLhwn81y4Y02twJ4F4BQbwMlBbdrfTBQMaVEC0Vo72nSPoY67tL+ZOKeesKVXI2VHG9ne/PpT9VqNnomRzc/JT516/OBb/ExIjn46DpdFM78QmSLlFyJTpPxCZIqUX4hMkfILkSlSfiEyRcovRKZI+YXIFAX5rAAXolhqrV6oAigOyOBrpDrncLCNzYfHVBe4yER8H24mxMUunJNYugeRHBRsROcsJYJ8ipJQmhQoU00Iz8+Qx1zsAoiDqfiqHDiUum8RlVZxYg8WwmdGxYqj9zZFm57zQ5d+sJR8/aCZX4hMkfILkSlSfiEyRTY/gGf/y99H22witM+4kGa7tvbPzTbZ0ZzE071RaEtXyY6sUlGJZK1KJ7uZj0kV6mC7mezVDtmrqU7AXKCyyAcS9+iNz6lUuAhmbPM32aavFNy3RFKWkY1vCWErVEjFOIGK5lZvJ4qobOL0q5lfiEyR8guRKVJ+ITJFNj8AJNa9jdaT283qquMUbGuz3exL8WcvF+SsNGk/2Z6RPQ+gqLZmqgAIF6/AEq9ZU5fe1H0LbPwypUa5gy432ODiF0BxN1z2T7QT/gqnQhwVeq8qKZuftkXr+OQ4Sa7zs+gbOB1r5hciU6T8QmRKofKb2S4z+66ZHTCze83s8t72HWZ2q5kd7P0+efjiCiEGRRmbvwXg/e7+YzPbBmC/md0K4J0AbnP3T5jZVQCuAnDl8EQdIgmbP4ptT8R2xyfRMUUx5Il1/sjmp2aQFYqP5/V4IF5/9grHyyfi1CskC9urNC6y77v3KTggYb+zbOzTSN83tNc7XASVbe2Er6UyF16jNkfPPdG0I/IDRG8/bUj1eSGb/+cffW/ioOFQOPO7+yPu/uPe308COADgDABvBHBD77AbALxpWEIKIQbPmmx+M9sN4MUA7gDwTHd/BOh+QAA4bdDCCSGGR2nlN7OtAL4O4Ap3P7aG8y4zs31mtm92drYfGYUQQ6CU8ptZHV3F/6K7f6O3+VEzO723/3QAh1Pnuvt17r7H3ffMzMwMQmYhxAAodPhZ1wNzPYAD7v7pZbtuBrAXwCd6v28aioRD4NnXUyJPKbcnJa2UCdjoo9lK5OArCPJBIqCHHW2VCsuaEmz1jj38jDx5Y3IS8n6eahLOu8JzUnCnYnK0RZ11luL3rjq/uoOvuhTfNhX4E+xPnMMMsSFPIWX+7c8H8GcAfmpmd/W2/RW6Sv8VM7sUwK8AvHk4IgohhkGh8rv797ByZOYrByuOEGKjUISfEJmixJ4V4OAZthPZNu+eVFxcswi2I7k7bGRnrq1B7fGrJrZxUA91mOXXknQbFET1cBxRXw8osYmKnEaFOPi9TOQBcWHUKgXk8Lh7HxpzYBT5WlKFO5JdljcIzfxCZIqUX4hMkfILkSmy+VeC7XeyI5NNHPqyv0M4KYfX+XmcIqrPyeNk5cxwGNeYKIgDSN6YdleLE4yKClqmbOTa0wXr+nRO6hq8Jl+bD8fVhUQSEidZUehDmeKcd//TxiXyMJr5hcgUKb8QmSLlFyJTpPxCZEoWDr9nf+4fwg11OiDlqCKnEDv4ynRwieCurancmChwhOXgi6RuVKZObtE51CmHD4+ShYoDVqKKOn1MPannHiXlFARKJR1+dE5tjjoHpRytXFm4EY75/W1tSVxjE9HML0SmSPmFyBQpvxCZIuUXIlOk/EJkipRfiEyR8guRKVms80fwOm9qWbwgOSRVzCMqtMFr8Ks3cU3eJyrmUSKxx6LKGyWKSnC3GdrAzXW8VpzoUtTBqEziS/Q8EglVRYU3ODYilYDF21IFO4vo1EPZmtPh/vbUJlbrTKCZX4hMkfILkSlSfiEy5YS0+c+69lPB2CbI5i1TjDIq5hHuTtmNhfZ4ZDfHh0Rx6HTfUsU8KKa8Uw9vXE0VHy26JsueaFoSPZOoYCc908R9omKb/DwSz702F46rSxSXX+K9i/0kxbSmwlewtI32T4fXvP+q9635HsNEM78QmSLlFyJTpPxCZIqUX4hMOSEdfhwY4lyIg45PNp9pFgT5pBx+7Fgq6OLqpYJN2HlV7JhqU6GN2ElY3KWXX2+nTGBUQVBT4fFIBDWVcfgtrF54g59Z6r2LRKNpsTUZv5g2OZJb28L7NE8aQDnnIaKZX4hMkfILkSlSfiEy5cS0+dnG56wUDlCpFgesRF1gEvZ8UUAK256dRKBMpb16gEoZmz8hWeERnIPDs0JU9LK4SXFhwE45m784Kac2v7pfpB8bv90ge34qPqe5NRy3ttBDmpLNL4QYQaT8QmRKofKb2aSZ/dDM7jaze83sw73tZ5nZHWZ20My+bGaNomsJIUaHMjb/IoAL3f0pM6sD+J6ZfQvA+wBc4+43mtnnAFwK4LNDlDXJcz56TbTNqClHba6ge2yquEVRZ9cSuTGR7clxAIlkkiieoMVJOcU2P9va8QGp7rjcpIP9InSJZJfecMj2Ou9PJkcV+E1Sz6y2yOv8qxcN8UTDkTYV4mjRGj6v6QNAh5t0NML7PrT3quicUaJw5vcuT/WG9d6PA7gQwNd6228A8KahSCiEGAqlbH4zq5rZXQAOA7gVwC8BHHX343PZIQBnrHDuZWa2z8z2zc7ODkJmIcQAKKX87t529xcB2AngPADnpg5b4dzr3H2Pu++ZmZnpX1IhxEBZk7ff3Y8CuB3AywBsN7PjPoOdAB4erGhCiGFS6PAzsxkATXc/amZTAF4F4JMAvgvgEgA3AtgL4KZhCrqifAmnUZUDcgq6tKYcflwNJw7YScgSBeSsfo1kt9jm6g4+3p/+vrU6qe7AhelOZSoZFQTxlAnYiRx+5OBLPbPqQriRq+iygy+VpMNVeVrTNE502G1PjFY13rVSxtt/OoAbzKyK7jeFr7j7LWb2MwA3mtnfAbgTwPVDlFMIMWAKld/dfwLgxYntD6Br/wshxhBF+AmRKWOf2FOmswp3cClTdIKTVJi0vbp6UQm2X6uJSrxFiTvGdnPq8IJnkshjKvZpcIBOyl+xVns9IXv0+kuY1e2JcA7r0H91azLcn0rSaW1ZvdsOB/R070sbquPlA9DML0SmSPmFyBQpvxCZMnY2/wuuCBN5qok16wrZvNHaeBnILu6QoZyyzdn/UJRwkpKL7WReo64uhYZzMjmoWqFjOO4h/szv1KNNq1ImKSn2E/A6f+IabPIXFAgBgNZU+E+wtI1sfFrXb0/G1+DCHB1K0mE/ApDoVCybXwgxDkj5hcgUKb8QmTJ2Nr8Qo8pD77xys0VYE5r5hcgUKb8QmSLlFyJTpPxCZMrYOfyqiyUKQhR1uSmRpBIV+IiiTxKyUeGNaExFJzhgpyvL6kE+lYXwBVs7cY1aKHylSYkv9Vh4rxVkMpFcUWXeBE5Vgs050SchR6Ky7nJak3FUV5S4Q0E9zW3h8akgnyioh4KeUoU7fCJVwnh80MwvRKZI+YXIFCm/EJky8jb/ee/4VDBuRF1c43OiTjkFRS9TtifDRSFTfgJOVKmQTV8le726mBDeWTi67zxVAEl9fFfCjc4+gMVEx55qwTzACTat+AGs9Rop+56fKxfq8EQlEi6q0eExFeJI+Tw61dWP8XriDa/J5hdCjCFSfiEyRcovRKaMvM0/+TivaxcXhKg2yRaLikwUXyOijL3KNj/ZxbbYpHGi80eH2/RSYY5meE5yXZz9BDUyaCvxZz6vyRdOC6lrFBQE6dS4yEjC9i6w8TuNlL+CrkHjqBhrqmArn8OFORLP46F3jVciD6OZX4hMkfILkSlSfiEyZeRt/saRsCpmpUk+gMR6M9vNkU3PdnWJdX6wbZ2weaPrtmgdn30NC2U6jtB96JrJSHi28dlPwPsBWKoS6nLo9ftE8bzRqbP9zvEHiTX7RnhMm+L22wmbv0PXSTciFYxmfiEyRcovRKZI+YXIFCm/EJkycg6/i8+9OhjX2uzgIycaj4E4OYYdcXx4wX4AsBo9Kg6KScGONZarDPT6vJUIDCIiyUiOZBENzqjh5CByCHYasVetPR1G+UQOP+6CVE8ECtFlizrwAnHhlajjMHdDTrx3TsVaogSiE9CJqJlfiEwprfxmVjWzO83slt74LDO7w8wOmtmXzSzRwVwIMaqsZea/HMCBZeNPArjG3c8BcATApYMUTAgxXErZ/Ga2E8BrAXwMwPvMzABcCOBPeofcAOBvAXx23RIdPRaOycZ18gF4soBlwcvytRdhiGztSiJQZmL1Lz8+yVUmEnJw8BAHLFHATtLn0Qht76jIRiPOwOk0wmfmZI+zjd+ajp9xYUAO2dpRkVTESTnswEgF+ayVZCEW2mZtkrVM8teYUXbmvxbAB/G7/LhTABx19+P/iYcAnDFg2YQQQ6RQ+c3sdQAOu/v+5ZsThyY/Gs3sMjPbZ2b7Zmdn+xRTCDFoysz85wN4g5k9COBGdL/uXwtgu5kd/+63E8DDqZPd/Tp33+Pue2ZmZgYgshBiEBTa/O5+NYCrAcDMLgDwAXd/u5l9FcAl6H4g7AVw0yAEas8+FoytTrYo2byVhP3KX0HiNXr6zOvDB2Bb4s4PTrJ0nrElHNdDg9YbiXVuSqAxKkZae2I+PCHVcITu05kKX397MmGvFxXRoAKm3BgDAJpbuPBGuD9qsJKqQ8K5UFHzlPgcbrDh/PJKuAnYxudksAff84Hii4wZ61nnvxJd59/96PoArh+MSEKIjWBNEX7ufjuA23t/PwDgvMGLJITYCBThJ0SmSPmFyJRNTey5qPLmaJtNhIEw7KyrbAmdaEgE9ETBNkXVa1PBNhw8Q9fw6anolPZ0eN/WVhpPkyOunkgwIdG4K3F8Qry/PclJOKtXxwGA5hQH4KxeHYe75ABAc3r1c9iZx42PU8cwqcCgqGMP+4CjhKJUl2Ian4CJPIxmfiEyRcovRKZI+YXIFCm/EJki5RciU6T8QmSKlF+ITNnUdf7qKTuibbZtazB2WrN3TvSZShSm4MKRPKY17MpCvLgcddilTkHtxH25wEXzJBpT4kuZwhS8zm9tXpCOz+FOtpy0w3IAQIvCFjgGIepimxC9TblOhWvlJYpqxAfEm1pbwofART4LO/Am+N8r3l94zLijmV+ITJHyC5EpUn4hMmVDbf6DP/k1Xr3r8t9tOPkZ0TGtHaHN35kIRWxtWT1uHYjt1aiLa4n4+epiaHzW5kObvzUZG7Rsay9Nc0HL8Pj2RKqaRThkWZtLxZ/X7Evg58Ex+F1ZwjEX4mA5kjb/FDW+KBI1sb8wtj9l80/TfevcgIRvcuIV4+wHzfxCZIqUX4hMkfILkSlSfiEyZUMdfue8YBe+te8zazrn/Ev+MRhz0Ql25nW3heN2otLsclIOv9pCeE59jj1e8XXYgRc5yTj4pETH2SggqVn8ec1BLVzcgp15qW3sAOTAmFS33E6DO93SNXj/oBxvU6GXsFILHX4Vkt0qcSRRp5PfPJjfKxZCAJDyC5EtUn4hMmVTE3vK8P2vjU+nlPPe8algHNna3AUnZfNzrVFyV3CAjiUKeBYW30zZ/JFfwFffnyiC2ZkgW5qOsTrb4n10SkpMV/VG2MWpVgt9AI0adTZOcOdrP7ZmWcYdzfxCZIqUX4hMkfILkSkjb/OPEz/81/EoAHHuh66JtnF8QYdjBWrkA9iSyMBphDZ8lcdki1cTNr/R2n+lUhwLMD2xFIpRDe+zpR7uF1008wuRKVJ+ITJFyi9Epsjmz5ADH3vvZotQmhf+518H40oiLr/Ixp+shuv8k7XmgKQbbzTzC5EpUn4hMqXU134zexDAkwDaAFruvsfMdgD4MoDdAB4E8BZ3PzIcMYUQg2YtM/8fufuL3H1Pb3wVgNvc/RwAt/XGQogxYT0OvzcCuKD39w0Abgdw5TrlESLg7td/dLNFOGEpO/M7gP8ys/1mdllv2zPd/REA6P0+bRgCCiGGQ9mZ/3x3f9jMTgNwq5n9vOwNeh8WlwHAmWee2YeIQohhUGrmd/eHe78PA/gPAOcBeNTMTgeA3u/DK5x7nbvvcfc9MzMzg5FaCLFuzBPFIIIDzKYBVNz9yd7ftwL4CIBXAnjM3T9hZlcB2OHuHyy41iyAhwCcCuA3g3gBG8C4yDoucgLjI+u4yAn8TtZnu3upWbaM8p+N7mwPdM2Ef3f3j5nZKQC+AuBMAL8C8GZ3f7zUTc32LVs1GGnGRdZxkRMYH1nHRU6gP1kLbX53fwDACxPbH0N39hdCjCGK8BMiUzZL+a/bpPv2w7jIOi5yAuMj67jICfQha6HNL4Q4MdHXfiEyZUOV38wuNrP7zOz+3vLgyGBmnzezw2Z2z7JtO8zsVjM72Pt98mbKeBwz22Vm3zWzA2Z2r5ld3ts+UvKa2aSZ/dDM7u7J+eHe9rPM7I6enF82s0Qngc3BzKpmdqeZ3dIbj6SsZvagmf3UzO4ys329bWt6/zdM+c2sCuCfAbwawPMAvM3MnrdR9y/BFwBcTNtGNXmpBeD97n4ugJcB+Mvesxw1eRcBXOjuLwTwIgAXm9nLAHwSwDU9OY8AuHQTZWQuB3Bg2XiUZV1fsp27b8gPgJcD+M6y8dUArt6o+5eUcTeAe5aN7wNweu/v0wHct9kyriD3TQAuGmV5AWwB8GMAL0U3GKWW+r/YZBl39pTmQgC3ALARlvVBAKfStjW9/xv5tf8MAL9eNj7U2zbKjHzykpntBvBiAHdgBOXtfY2+C93w71sB/BLAUXc/XltrlP4PrgXwQQDHa4WdgtGVdd3JdhtZw88S27TUsA7MbCuArwO4wt2PmaUe8ebi7m0ALzKz7ehGip6bOmxjpYoxs9cBOOzu+83sguObE4duuqw9+k62O85GzvyHAOxaNt4J4OENvH8/lEpe2gzMrI6u4n/R3b/R2zyy8rr7UXRrPrwMwHYzOz7xjMr/wfkA3tCrWnUjul/9r8VoygpfR7LdcTZS+X8E4Jye97QB4K0Abt7A+/fDzQD29v7ei65tvelYd4q/HsABd//0sl0jJa+ZzfRmfJjZFIBXoetM+y6AS3qHbbqcAODuV7v7Tnffje7/5n+7+9sxgrKa2bSZbTv+N4A/BnAP1vr+b7CT4jUAfoGu3fehzXaakGxfAvAIgCa631IuRdfmuw3Awd7vHZstZ0/WP0T36+dPANzV+3nNqMkL4AUA7uzJeQ+Av+ltPxvADwHcD+CrACY2+5mS3BcAuGVUZe3JdHfv597jurTW918RfkJkiiL8hMgUKb8QmSLlFyJTpPxCZIqUX4hMkfILkSlSfiEyRcovRKb8PxNDioo/YLA+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Mask zone of raster\n",
    "zone_raster = np.ma.masked_array(data_raster, np.logical_not(data_mask))\n",
    "plt.imshow(zone_raster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1185.0, 1368.0)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.ma.min(zone_raster), np.ma.max(zone_raster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 70, 146, 270, 403, 426, 391, 394, 314, 134,  53]),\n",
       " array([1182. , 1201.3, 1220.6, 1239.9, 1259.2, 1278.5, 1297.8, 1317.1,\n",
       "        1336.4, 1355.7, 1375. ]))"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.histogram(zone_raster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1268.6926229508197,\n",
       " 1268.6926229508197,\n",
       " 1268.0,\n",
       " 38.405802119545115,\n",
       " 1475.0056364456561)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.ma.average(zone_raster), np.ma.mean(zone_raster), np.ma.median(zone_raster), np.ma.std(zone_raster), np.ma.var(zone_raster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 栅格图像分区统计(tif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fcfbb671e10>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQMAAAD8CAYAAABzYsGzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAADDRJREFUeJzt3V+o1/d9x/HnayfaNGuH0ZjgVGYGdqQXiYFDlpFdZLbZXFaaXmSQrAwHgjcbpKzQmA3GCrtIbprcjA5ZQr0oTbq1Q5GCE6uMwjA5JsZqbaIN6SpKtFLpRiFM+97F+ZidnRx3fuec3/f35/h8wOH3+37O9+v3jXzzyvv7+X0/P1NVSNKvDLsASaPBMJAEGAaSGsNAEmAYSGoMA0mAYSCpMQwkAUsMgyTbkryV5GySXf0qStLgZbFPICaZAN4GHgHOAa8BT1bVD250zB2rJ2rTxhWLOt9y8faJ2xZ8zCfu/UUHlehmcezE+z+tqrXz7XfLEs7xAHC2qt4BSPIy8BhwwzDYtHEFrx7YuIRTjr8/+PUtCz7mwIHjHVSim8XEurM/7mW/pdwmrAd+MmP7XBuTNIaWEgaZY+xD9xxJdiaZSjJ16fK1JZxOUpeWcptwDpjZ828Azs/eqap2A7sBJu+71SWSi7CYW4u5HDjv7YZubCmdwWvA5iR3J1kJPAHs609ZkgZt0Z1BVV1N8hfAAWACeKmqTvWtMkkDtZTbBKrqO8B3+lSLpCHyCURJwBI7A42XfkxEOgm5fNkZSAIMA0mNYSAJcM6gc/16YGhU+ADU8mVnIAkwDCQ1hoEkwDkDDYnPPIweOwNJgGEgqTEMJAGGgaTGCUSNLSch+8vOQBJgGEhqDANJgHMGfbfcFiYtdy68+l92BpIAw0BSYxhIAgwDSY0TiFIfLIcHoOwMJAGGgaTGMJAEOGcgjYxhPwBlZyAJMAwkNYaBJMA5gyVxUZJG0Yevy7M9HWdnIAkwDCQ184ZBkpeSXExycsbY6iQHk5xpr7d3W6akrvXSGXwN2DZrbBdwqKo2A4fatqQxNu8EYlX9W5JNs4YfAx5u7/cAR4Cn+1iXpEWa/dDRxLrejlvsnMFdVXUBoL3eucg/R9KI6HwCMcnOJFNJpi5dvtb16SQt0mLD4L0k6wDa68Ub7VhVu6tqsqom166ZWOTpJHVtsWGwD9je3m8H9vanHEnD0stHi98A/h34rSTnkuwAngUeSXIGeKRtSxpjvXya8OQNfvWpPtciaYh8AlES4EKlBXFhkpYzOwNJgGEgqTEMJAGGgaTGCURpjPXzX2GyM5AEGAaSGsNAEmAYSGoMA0mAYSCpMQwkAT5n8P9yYZJuJnYGkgDDQFJjGEgCDANJjROI0hjp58Kk2ewMJAGGgaTGMJAEGAaSGsNAEmAYSGoMA0mAzxl8wEVJutnZGUgCDANJjWEgCTAMJDWGgSTAMJDUGAaSgB7CIMnGJIeTnE5yKslTbXx1koNJzrTX27svV1JXenno6Crwxap6PcnHgWNJDgJ/BhyqqmeT7AJ2AU93V6p08+nyy0xmm7czqKoLVfV6e/+fwGlgPfAYsKfttgf4XFdFSureguYMkmwC7geOAndV1QWYDgzgzn4XJ2lweg6DJB8DvgV8oap+voDjdiaZSjJ16fK1xdQoaQB6CoMkK5gOgq9X1bfb8HtJ1rXfrwMuznVsVe2uqsmqmly7ZqIfNUvqwLwTiEkCvAicrqqvzPjVPmA78Gx73dtJhR1xlaL0f/XyacJDwJ8C309yfWrzr5gOgW8m2QH8B/DH3ZQoaRDmDYOq+h6QG/z6U/0tR9Kw+ASiJMAwkNQYBpIAw0BSYxhIAvx2ZGlkDHJR0lzsDCQBhoGkxjCQBBgGkpqbYgLRRUnS/OwMJAGGgaTGMJAEGAaSGsNAEmAYSGoMA0mAYSCpMQwkAYaBpMYwkAQYBpKaZblQyYVJGgfD/maj2ewMJAGGgaTGMJAEGAaSGsNAEmAYSGoMA0mAYSCpMQwkAYaBpGbeMEhya5JXk7yZ5FSSL7fxu5McTXImyStJVnZfrqSu9NIZvA9srar7gC3AtiQPAs8Bz1fVZuBnwI7uypTUtXkXKlVVAf/VNle0nwK2An/SxvcAfwt8tf8lzs+FSRoHo7Ywabae5gySTCQ5DlwEDgI/Aq5U1dW2yzlgfTclShqEnsKgqq5V1RZgA/AAcM9cu811bJKdSaaSTF26fG3xlUrq1II+TaiqK8AR4EFgVZLrtxkbgPM3OGZ3VU1W1eTaNRNLqVVSh3r5NGFtklXt/UeBTwOngcPA42237cDeroqU1L1evuloHbAnyQTT4fHNqtqf5AfAy0n+DngDeLHDOiV1rJdPE04A988x/g7T8weSlgGfQJQEGAaSGsNAEmAYSGoMA0mAYSCpGbt/UclFSVI37AwkAYaBpMYwkASM4ZyBNA5G/YtM5mJnIAkwDCQ1hoEkwDCQ1BgGkgDDQFJjGEgCDANJzcg/dOTCJGkw7AwkAYaBpMYwkASMwZyBNA7GcWHSbHYGkgDDQFJjGEgCDANJjWEgCTAMJDWGgSTAMJDUjNxDRy5MkobDzkASYBhIanoOgyQTSd5Isr9t353kaJIzSV5JsrK7MiV1bSFzBk8Bp4Ffa9vPAc9X1ctJ/gHYAXy1z/VJI2c5LEqaS0+dQZINwB8B/9i2A2wF/rntsgf4XBcFShqMXm8TXgC+BPyyba8BrlTV1bZ9Dlg/14FJdiaZSjJ16fK1JRUrqTvzhkGSzwAXq+rYzOE5dq25jq+q3VU1WVWTa9dMLLJMSV3rZc7gIeCzSR4FbmV6zuAFYFWSW1p3sAE4312Zkro2b2dQVc9U1Yaq2gQ8AXy3qj4PHAYeb7ttB/Z2VqWkzi3lOYOngb9McpbpOYQX+1OSpGFY0OPIVXUEONLevwM80P+SJA2DTyBKAoa8UMlFSdLosDOQBBgGkhrDQBJgGEhqDANJgGEgqTEMJAGGgaRm5L4dWRo1y/WbjWazM5AEGAaSGsNAEjDgOYO3T9zm4iRpRNkZSAIMA0mNYSAJMAwkNQOdQPzEvb/gwIGlPcDhBKTUDTsDSYBhIKkxDCQBY7hQadwWjTjHMX7G7RrrFzsDSYBhIKkxDCQBYzhnMG7G6f7T+Y2bm52BJMAwkNQYBpIAw0BS4wSiPjBOk53ghGe/2RlIAgwDSY1hIAmAVNXgTpZcAn4M3AH8dGAnXppxqhXGq95xqhXGq96Ztf5GVa2d74CBhsEHJ02mqmpy4CdehHGqFcar3nGqFcar3sXU6m2CJMAwkNQMKwx2D+m8izFOtcJ41TtOtcJ41bvgWocyZyBp9HibIAkYcBgk2ZbkrSRnk+wa5Ll7keSlJBeTnJwxtjrJwSRn2uvtw6zxuiQbkxxOcjrJqSRPtfFRrffWJK8mebPV++U2fneSo63eV5KsHHat1yWZSPJGkv1te5RrfTfJ95McTzLVxhZ0LQwsDJJMAH8P/CHwSeDJJJ8c1Pl79DVg26yxXcChqtoMHGrbo+Aq8MWqugd4EPjz9vc5qvW+D2ytqvuALcC2JA8CzwHPt3p/BuwYYo2zPQWcnrE9yrUC/F5VbZnxkeLCroWqGsgP8DvAgRnbzwDPDOr8C6hzE3ByxvZbwLr2fh3w1rBrvEHde4FHxqFe4DbgdeC3mX4w5pa5rpEh17ih/Qe0FdgPZFRrbfW8C9wxa2xB18IgbxPWAz+ZsX2ujY26u6rqAkB7vXPI9XxIkk3A/cBRRrje1nYfBy4CB4EfAVeq6mrbZZSuiReALwG/bNtrGN1aAQr41yTHkuxsYwu6Fga5hDlzjPlRxhIl+RjwLeALVfXzZK6/5tFQVdeALUlWAf8C3DPXboOt6sOSfAa4WFXHkjx8fXiOXYde6wwPVdX5JHcCB5P8cKF/wCA7g3PAxhnbG4DzAzz/Yr2XZB1Ae7045Ho+kGQF00Hw9ar6dhse2Xqvq6orwBGm5zpWJbn+P6VRuSYeAj6b5F3gZaZvFV5gNGsFoKrOt9eLTAftAyzwWhhkGLwGbG4zsiuBJ4B9Azz/Yu0Dtrf325m+Nx+6TLcALwKnq+orM341qvWubR0BST4KfJrpybnDwONtt5Got6qeqaoNVbWJ6ev0u1X1eUawVoAkv5rk49ffA78PnGSh18KAJzkeBd5m+l7xr4c96TJHfd8ALgD/zXQns4Ppe8VDwJn2unrYdbZaf5fpNvUEcLz9PDrC9d4LvNHqPQn8TRv/TeBV4CzwT8BHhl3rrLofBvaPcq2trjfbz6nr/20t9FrwCURJgE8gSmoMA0mAYSCpMQwkAYaBpMYwkAQYBpIaw0ASAP8DQzz8WkhIo+gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "raster_region = gdal.Open(\"0.tif\")\n",
    "region_band = raster_region.GetRasterBand(1)\n",
    "region_array = region_band.ReadAsArray(0, 0, region_band.XSize, region_band.YSize)\n",
    "\n",
    "plt.imshow(region_array)\n",
    "\n",
    "# region_band.GetNoDataValue()\n",
    "\n",
    "# print(f\"{region_array[region_array == 0]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "255.0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "region_band.GetNoDataValue()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "region_x_origin: 12489296.69225333, region_y_origin: 4556772.136311595\n"
     ]
    }
   ],
   "source": [
    "trans = raster_region.GetGeoTransform()\n",
    "\n",
    "region_x_origin = trans[0]\n",
    "region_y_origin = trans[3]\n",
    "region_pixel_width = trans[1]\n",
    "region_pixel_height = trans[5]\n",
    "\n",
    "print(f\"region_x_origin: {region_x_origin}, region_y_origin: {region_y_origin}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "51, 50\n"
     ]
    }
   ],
   "source": [
    "region_x_offset = (region_x_origin - x_origin) / region_pixel_width\n",
    "region_y_offset = (y_origin - region_y_origin) / region_pixel_width\n",
    "\n",
    "print(f\"{region_band.XSize}, {region_band.YSize}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1278 1295 1300 ... 1348 1358 1374]\n",
      " [1276 1293 1301 ... 1347 1359 1375]\n",
      " [1275 1291 1299 ... 1348 1360 1371]\n",
      " ...\n",
      " [1218 1233 1240 ... 1349 1354 1358]\n",
      " [1221 1236 1244 ... 1335 1340 1346]\n",
      " [1225 1238 1246 ... 1323 1328 1333]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fcfbb4b66a0>"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQMAAAD8CAYAAABzYsGzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnVuMJdd1nv9VVefS3XMhhzOkGA4TygBjyA8xBUwUBcqDQlsBoxiWHxTAshHQAAG+xICUOLDIBAhiIA/Si6WXwAERCeaDYcqOHZAgDDgEIyIIEFAaSZRMmpJJKZI1GJpDihpyZrr7XKp2HrrIPvtfa7qqz8yc7kn+DyCau7ouuy6zutZf62IpJQghRHHQExBCHA5kDIQQAGQMhBAtMgZCCAAyBkKIFhkDIQQAGQMhRIuMgRACwDUaAzN7wMy+Z2avmtkj12tSQojVY8tGIJpZCeCvAHwMwDkAXwfwqZTSX15tm3JjI1UnTuwuGPhjm+09n3C6yfJxTeNwGzpuj9060xnM1cp8GZ9PVTZum9LyZRWNh8U82Cbfb+Mm66npBJa59Ty3nWU1zS1fpwhuQIF8HT6fCtFx8nPkvU6D89lOg2w8S2U2ntN4Z7/5cSx8gPaG9wH4692443j4Wo5tlo1HRT4GgA3a0Te/M30zpXTq6rPdoepaYQ8+BODVlNIPAMDMngDwCQBXNQbViRO461995r1xc/vUrzPIH6yC/vHUc3/z6iktu5SfVjHxl9nmdCP8s4dEV6dez1dKI7/R4Eh+Tmtr+fjExqbb5vhwOxufGl/OxqfHP3XbHC3zbbab/KHnBw0ALtXjbDxpum9/Sf8obxtcceucHFzKxreU+Tr8AAPARjHJtym2aOyfjZNlfp9nKZ/b+do/G9+d3pGvM7s1G785O+q2YYPB16Do+IMV7QMA5k1ujLeaYTZmowoAx6r8Pv/d8Wv5ePi62+bvj/J7P/xbP/jR3rPd4VrchLsA/HhhfK5dJoS4CbkWYxC+WbuVzB42s7Nmdra+4v+qCCEOB9fiJpwDcPfC+DSA87xSSukxAI8BwOj03ckW/PnEvj0AG+39CjYc+VdOfqGcr+X7rStv84pJvix4Q0PizXhqwVTns/z1cFLml/idcuQ3IoZlrhHcMQpcI7LFY/Id+fdA4Lv73ToGpFecqi65dd5XvZ2Nbym8K8Ss03xLmlv0V+qtOr9JbzT5tfzx7Da3zRvzY9mY3YI3pkfcNtMO9ynSTZi5e3iASZ3vl125SFthN6GhK1MEcylt4Jb14VreDL4O4F4ze7+ZDQH8KoCnrmF/QogDZOk3g5TS3Mx+E8CfY+dvzJdTSi9dt5kJIVbKtbgJSCn9GYA/u05zEUIcIIpAFEIAuMY3g6VY1EgCAbGpc/vEgmEZBO2Mxvk6LMM0HIcAoKFDWzAX3lGqekQq8XHofObBt/BEAUNTEppYHAQCwY2FpEDAKnuIjMw6fe8/Vb3j1jlV5qLiuvkgqS6GNP9LHOQBYEbn9N3Jndn49flxt83fTPJlLBi+M11z20yb/B5VRT43jheIiILA+D6zyHhqLY8vifbjgrUC0fFys+2W9UFvBkIIADIGQogWGQMhBIAVawYGygFgxx2Bv0/+2XoQdNQkSpQhH2868Kc5rfJlae7tYnLCAvnpldcvOFGppHWq0kc38XzXKx+X38VmnQfgRLHxXXC+AxDkJhQ+ipTzCI4W+XXbDjKialr0FsXpX0n5GPBBRd++8rezcRRAxJrApVl+nbbn/tnoSvqqA82A9YAmCEijy4KNYX7dpoGexHkQnNAV6T6bKYig64HeDIQQAGQMhBAtMgZCCAAr1gwS6PN3ERS9IJ9oMMi/WQ8r/w17SH44+17NyPvgm4PcJ51MA9+x3ttWcq0FwMdBjIe5xnFsnOfxA8DxYZ7Lf3KUf2+OElg2ycfmegZb9XLJKl1wsZCdZfn1LZqafu+v4yXaz1t17u9fmPs6A+emuWbw2nYeQ7A59zrDxQlpBhO67zN/Pixx1B3PAeDvuwWywzo9h1z45sggeDaqvZO+Nhuf+Hap8fvpg94MhBAAZAyEEC0yBkIIADIGQoiW1SYqFUC9tiua2NALcCzKjapcjDoy8GIgB+3cOspFlyYQsN6u8uKgVwZefKqDoKj8uF7YYzHz2CgP5Dk29IE97xvnyT8jqi4UBRDxsnfm+fls1f58Ls33rrL0RlCF+QgFQJ2ofDINi4pcADUKjLlCwtfFej0bs6AIeNGUE33enubXAADe3qLrMqFqyYFw7ILNmOD31Si/dhvrXsTjZ4Of01NDf22ZmoKbfhJcp/VCAqIQ4hqQMRBCAJAxEEK0rFgzSEjru37TYOATKkbD3PcaUZDRWuUTlY5R0M5ama/TVe0WAEal95c5+YRhrQKIA0cWOT7Ycss4GYjH7CsDwCXSCC5SQk4UgLM1z/1l9rkHhb8fl0kz+NHwpFuH/X9ukDINCpVwMtClhhu8+GAgDqTapnGU6DOjZXNqwtPMgr+H0bJFBkGwmUtIC56NYX5duHlOH1hL4YAvwDex6YveDIQQAGQMhBAtMgZCCAAyBkKIFhkDIQQAGQMhRIuMgRACgIyBEKJltdWRy4TBkd0giSiZY4OqwXCgRlQ5mJNphpRwwxVlAaAa50khk2b/1YRHQZBOF1HSVE02+eI8T9q5HAQQXaYqvxxktD33QTuTICgnJ9omf0S+v+mDjn46zOfLiVaDoN+96xRE9yhKzuLrxBWgejS48gRVsV2nL95x0H2LKx1xsBzgn11u7c6dmwDgbXoWuLtWdG2jalR90JuBEAKAjIEQokXGQAgBYNWagSUMF4pADCrv7wyoAMSYEoiqwE/ngg/sj1aBX1VT+doTw/0XhOiTANWHt2d5khH76ZyEBPikI+4MNJ13J+0w3L0KALbpfkQJXZwwtNGh4QBeV2DNIOpqNKHrzd2SgyZGrtJxwxpB4P8bFS9J1CUr+hNaFN3+P+tFnCjGyWeAv3aspUSVs7crrzH1QW8GQggAMgZCiJZOY2BmXzazC2b24sKyE2b2jJm90v689cZOUwhxo+nzZvD7AB6gZY8AeDaldC+AZ9uxEOImplMBSyn9TzO7hxZ/AsBH2/9/HMBzAD7b54CL1YGi6sMshA2LXJxqgtCSbapsxMEcESyERYFJR6v9V6LpamvGYiHgA4a4Ys/bEy8sbU+pyq+r6NNt57vaxwFAXZHINd27wjIAzEmI5KrAgBfGOICLxTXAC8XXhUh1pMfH6JlLgejI4ussEGNZ8ORqzpPS/3NkwZyreB2GoKM7UkqvAUD78/Yl9yOEOCTccAHRzB42s7NmdrZ+Z+8mkkKIg2NZY/C6md0JAO3PC1dbMaX0WErpTErpTHls/WqrCSEOmGWjZp4C8CCAz7U/n+y74XzBT60Dv4p9rzn5wlFL9gEFeHDV4qjqL/uxUQJUFCyTzTXwYTmw5CeT3ABemXmfm4OMNqlN+JVtH0TCnYD6tA3v7BTUg+ncPzLbdL25qvR27TWcId2TYbl3EBLgr3d0/a8LfJ24c1agGfD1j67TZpHfR54/B28BPnhpXuXBcawhAL5adV/6fFr8QwD/G8DPmtk5M3sIO0bgY2b2CoCPtWMhxE1Mn68Jn7rKr37hOs9FCHGAKAJRCAFgxYlKTWOYLnTBbaI4A/LP2HfkLs2ALyzhxoH/yUlSG0OvGXASCCfgRN/Cr1DMAMcQvBN8p+eYAR5Pt4OORBxHwP5z4Nd2YSOvrbD/Pw+0CY4rMKNkmuA6sb/M3+AjouSf/WL0LERhBo55Ptfoys5n+fy3yu5v/VxsJtJJWA/jax0ljr05O9p57Ai9GQghAMgYCCFaZAyEEABkDIQQLattyd4Y5lcWhJUeIlfNVWYCkYVNmlF7bOOgEQADav0+CyoDcft3FgOLICGK17lMguHmxAcQbU+oVfqU2oZvBbeJxFdbQjBMQ5p/IOj2Edg4sYq3iWZWFyw6dh+pqxpSryAkXqXPJjy1qd9oPsmvQXQ+HJjkhO5AIK0oOM4lgQUBdRcGEhCFENeAjIEQAoCMgRCiZeWaQXFl17eyvfOA2pXyYRC/4peRRtAMvP82GeY+Hif+AP0CkxguMrI1y/fL+gAAzCb5Oon9z4k/aZsvoRFUeyfcpB5PQ1hghPQWLlpTBJtMyaeO/GWGNQPX6ChIfKs6AtKaQE8C6VS2TfpMoK00dM98+hBQs5ZV7h1gB/ggO07ciwKxzpfHg6N3ozcDIQQAGQMhRIuMgRACwKo7KjVAubnrbxWzqJvN3vtw3W0AcMPeRN2S6vWoUAaNg2NxTAAn7UTwfidUqIQTWgAgTckmkx5QbkZCST5kdzNoYuzCCJyGEJwfF/+czfwjU5Pf6hKVAr+Wj+QS0oJtWHvgb/B9cIlukWZQ7x3DET2jxVZ+j/g5AICm4oCYfCUXUwOvK/QpYrssejMQQgCQMRBCtMgYCCEAyBgIIVpkDIQQAGQMhBAtMgZCCAAyBkKIltUGHdXA4PJC0FGQ99MZdFQFgTEcs0NBRikoRpEGlKQTHHdGAUKJA3165Aq5fURdjZq9g4zKvIlOuyOafxQ8Q3DDXk64iboL8/Wug+7OUYBTTpBoxU2LOgKXAJ/M5ALHgkQll9y0RBMmd92CWKfSJTMFFaGH/ADlwzCgjtaZjboLukQJW33Qm4EQAoCMgRCiRcZACAFg5cVNgGpzd1hOgsQM8t1ZD2gCzWC+cR3mFvherjDpEkki3DXKJSUBMFpWTGzP8c7CfJgoiydKenGVSlkniTYhjaAOtAn2W5fpjWxGBVKD43ARVfaNo+QmJurixXDzbTcOEuyaHs2eOOHJ6zyRZkPjebeGsxks64PeDIQQAGQMhBAtMgZCCAAyBkKIltUGHSWg2toVTaKgI1fpmAWUqIIPLWtoGxfsgSDAIxDcGgoYaqKStwyXHOJAnkDcMRKkWKAKmua4IKnkGzX547B4RvuwWRBQRJE9TVDBtzOQp09AFO832CeLipF4xnAXo4bO0Sb+gXKVjbgle1DVm4suhRWUqbKUceBYdB1ZQC9pm+1AdCy729tH6M1ACAFAxkAI0dJpDMzsbjP7qpm9bGYvmdmn2+UnzOwZM3ul/XnrjZ+uEOJG0UczmAP4rZTSN83sKIBvmNkzAH4DwLMppc+Z2SMAHgHw2b12ZA1Qbu+Oizrw08knmlMDoqijkgvA4WrJXAUY8H5s1MW4q7Nx5Ap3bML6AAAUlEjFWkrYearce52iR3XkgjSQsPozF1AuAl2BF/SJOiKNwN3XHtqEi5mKukxxt2rSDMor/nzKLUoco/vB1w0I/P3okaP77Cpa97huHIQUBeH1SVqL6HwzSCm9llL6Zvv/lwC8DOAuAJ8A8Hi72uMAfmWpGQghDgX70gzM7B4AHwTwPIA7UkqvATsGA8Dt13tyQojV0dsYmNkRAH8C4DMppXf2sd3DZnbWzM7Ot68sM0chxAroZQzMbIAdQ/AHKaU/bRe/bmZ3tr+/E8CFaNuU0mMppTMppTPV+HpkFAkhbgSdAqKZGYAvAXg5pfS7C796CsCDAD7X/nyyc18JqIJMxUU4YIjp05KdBcOogkwffLAJl6YJtumRvcZ0ZskFAiJXXXKiXSBU+gl3V3vybcKCG9Cls0aC1hIlh9zse7Q943VYxFsUtN9bRpWlWNCNBESmzzp9npWGBXQObgqejWa4TN5ov68JHwHwLwD8hZm90C77t9gxAn9kZg8B+GsA/3ypGQghDgWdxiCl9L9wddv/C9d3OkKIg0IRiEIIAKuudJSAYrbr+TWDbt8mcWJGEEzDmkB3td54bgxXGOKAoSAuJg5E6sAFHfXRDPbfjRz+BY+SkMLYrG5dIQwEWzxqoA/0CbAJJpMNS2qDHlUtdq3raR2uagx4zcBpCPPuCl2cPLez0t5ziyp/O12B5h9d+ug+9kFvBkIIADIGQogWGQMhBIBVFzdpUlYROUqyYF/LxRAE5otjE5byR4NiFBxXwN+OowIWy8Dfiv048FFpzNWRo8/4vIx9+SIQPJxEEFzcrryYUFPgkI0+BVCoGnJ1pTvOwBVwoTEnIQE+9qDczucWFZtxGkEP3Yo1s1AP4+efk7WCbeq1G5SoJIT4/wMZAyEEABkDIUSLjIEQAsCqg44A2EKGTR9hySVqBDPuCnoJ6SH+cYBK0aNKrtsHB5YEgg8LUixyRUFHXuvrbtXVRRNsw6JiNP8w2KdjGxbCrMd8+Tq4YKAgOcjNjS5TdG2rTQpumtE1iNRZfk5Lv858jRb0EMc7BcSRFwvn6xIQhRDXgIyBEAKAjIEQokXGQAgBQMZACNEiYyCEACBjIIRoWXlxk8Xvvr2Sjnp8i+0yaVFCERdE4cKZO8voMJyoFHxf7+yS06OIqoszWKq4Zo8ipK6Lkd+EuyxFhTP4OrhrEMVJdIUVRIVWOM6AE4qCYrsujoCvbXAPi9ne3+nNVaMF6mF+MedBIfD5mKZCMQLRs+26XnGcQdB9O230CIAJ0JuBEAKAjIEQokXGQAgBQMZACNGy8kSlRbg6D9BdDTnsjsQVcFnA6lPFOKiS4xKTlhD2Cq6aG1xxn5iU9hxH8LUME7p6zJdhwbCMRFOav0tC6pMzxfcwvB/5uNqi6xRULeJtjE4oEhC7mK35E5odzcdRslBNAmIdJBkxTgClTaJEJRsu0dYLejMQQrTIGAghAMgYCCFaDqCj0q4/E3Y+YnfMBe2EkTH5Ji7mItqG1uAOy+E6+TjSDLr8+zC4iTUD7vrbI4akGfA1CDoFdRVAiQKKelSe7gy0CuBtXOBSD22Cg46ia++6U9U99Bea/3yc/82sR36b2Ua+39kxf5y0Ric16Pbt6+38H4kLoItOZ8mq3XozEEIAkDEQQrTIGAghAKxaM7C8i1LkW7oOMT1qO/K34oY+bEf+J+sVvRKV3Ddrv1+XdNTxnbvPcfrEGXS2+AVQ0zqF62IU3BCaWyjZ0J8Ujq2Ir1Pac51QM6BlFXc6ijSDDo0gKm5aj43G+e/nG8E26/kJOH0AgI3zZWXVrRnMu7pgBzfEenSnitCbgRACgIyBEKKl0xiY2djMvmZm3zazl8zsd9rl7zez583sFTP7ipkFmdVCiJuFPm8GEwD3p5R+HsB9AB4wsw8D+DyAL6SU7gXwUwAP3bhpCiFuNJ0CYkopAbjcDgftfwnA/QB+rV3+OID/AOD3Oo+4INiE3XlYIOHiPEEADsdYVNPuoIteAiIn4PQwnS4wicZOAAqOw0JYVMGnax+RGMhakxP+AtGR1+FgoWhZVxcjIOgi5bpKRWJgPq4m+UlHFYqcQMgBReuBGDiyvcfBO7BLGApEvIIEw9G4O3OsGuQn3dDDPt2kVk7AjQ06MrPSzF4AcAHAMwC+D+BiSu8VtDoH4K6lZiCEOBT0MgYppTqldB+A0wA+BOAD0WrRtmb2sJmdNbOzs+mV5WcqhLih7OtrQkrpIoDnAHwYwC1m70WunwZw/irbPJZSOpNSOjMYBlUihRCHgk7NwMxOAZillC6a2RqAX8SOePhVAJ8E8ASABwE82Xk0A5rB3poBd112yUJBMEpJ/n5XtWEgKAYSmMWuYibLJMakwJ/jLr/s+4ZBR65iS7dd9+dIvnCwDc82KgbiqyNz5elgG1eluMc2tN+Sgo4WA9p2N8qHs3WqYjz228woqGhGf8OiisTusFXwbJC4MqBKMYPKn3RNz8u8zv/RNLW/78uFHPWLQLwTwONmVmLnifujlNLTZvaXAJ4ws/8I4FsAvrTkHIQQh4A+XxO+A+CDwfIfYEc/EEL8P4AiEIUQAGQMhBAtK81aTIVlAk4k+LgMPg6U6VHNtpzQPiP9rSMAJ6JPZp2rbOwyEr2a2VkdKdjGnxOv40+o5AAiF+wUHJuPEwYQ7S3+Re3IOoOMehT4XRSjAaAe+ueJBUIe11Gl4yN8nHxcj4PgJm7X59YAijI/KRYMh0Hp6dEovyl1k99EFiUBYHN7ucwAvRkIIQDIGAghWmQMhBAAVq0ZWO6bx74yZ5LQPsK21Ry0s8Tkogo+nMzELcGDDj4VJRVxkpEF59xVXbicBQ600y+4ulPkp+cXL0VBOl2EiUp7+/thchPpIH0qLPM6k/X88Z0e8Q+Hq1LEGkIQFMsaQTOkwLEeVb0tqGI0HOYP0Pogf1CPjajcM4Bxma8zrfNztuDizupogt3ozUAIAUDGQAjRImMghAAgYyCEaJExEEIAkDEQQrTIGAghAMgYCCFaVtteDXnQkG+dDpQclMPBHGFyEI8p0CeKberRWYwDX3g/UTVeDjIqt6mlVo8AIjbRxZa/UDyXpqIElpk/oWbOwTM9go74nOvuDCKuSBwlKvXIq3LMxyWN926DBgCzI5yoRNOIKh1Te3sOOnKVkOFbp3FVYwAYUWLSuMoDitYrH8V2tMqz7iZlfg2mrh+hr6DUF70ZCCEAyBgIIVpkDIQQAFasGVjKE3miqsVum44qukBQUIR9+aiLUY8Sslx8pU/V32KaH6zcntPvg434nMiXL7Z85hX75SVtkwpv5wvSFaJkIHccqogSFWfp2g/vA4DXaKgKcD3yvjBXtOZiJlFxE+5+lCoeB4VKaB1OXEoDv001zO/raOTv2cYw1wSODfPEJNYHAGCDlg2afHLbA79NVGW5D3ozEEIAkDEQQrTIGAghAKxaM2gSqq1dn7MKWi+6jr5cBCPw/4sZr8OVMoLvwoHvyzhfl+MOokKldGwjzcCmQeWVhoWR/Lg2C+IMjDUC0jfKoMAFV5ctWEPoISJEsQmDrmIagc5A1zYN+HyCvVA3ZC4y0qfoCO83Og4XNwWPgw7LJccQDPw947gCLlzC+gAAHKHqvjMKqpkFcQasTfRFbwZCCAAyBkKIFhkDIQQAGQMhRMtqBcQ6YXhxQTRhoQ+AUSKMOXEwEP5cRx8KugiOE4mKfjIsuNE42i8d24mZkYDYddx5kNHFm1R0K6P2SCQqpiofh/IhC5XlIForo2FBka8bfGXmepD/XarHXhirSWRsOIBouaLA+ye4UCV3SwqShbhj0ojabbFYCADHq81sPKEIqAnfdwDHh77Kch/0ZiCEACBjIIRokTEQQgBYtWYwb1C9tRBpFBTKsDn5Wuz/R8U1eJ0eekDqsY5xsk8UyMNwoFKkK3RBGkGKdBLCWFcIfEkXLNMR7LSzH9IZhn6/NXU24kIrUSJTQxoBF1qpx34j17V7icI3bi5RURsXwEW/X6YTFYCCKsUMSDMYBCdwtMj9//WiO6Bour7cP2u9GQghAMgYCCFaehsDMyvN7Ftm9nQ7fr+ZPW9mr5jZV8wsqCYnhLhZ2I9z8WkALwM41o4/D+ALKaUnzOw/A3gIwO/tuYf5HHjjrd1x4IMn9n0pSSfy9d03doY1hZ6kAe2XC32OgwqcvI+1Ub4g8v/ZV2841iKIGWDthOYaFhThdUb5N+sUJByxbz8P/FGOCahHe/v2gPfdWTOIdIaoeMl+4UQ37l69s05HEZs+ISrBsoIOXtGO14M4gzG1FGcNYaPw24yWakPe883AzE4D+GcA/ks7NgD3A/iv7SqPA/iVpWYghDgU9HUTvgjgt7Gbi3obgIsppXf/ZJ0DcFe0oZk9bGZnzezstFkuMkoIcePpNAZm9ksALqSUvrG4OFg1fHlKKT2WUjqTUjozLLpfq4UQB0MfzeAjAH7ZzD4OYIwdzeCLAG4xs6p9OzgN4PyNm6YQ4kbTaQxSSo8CeBQAzOyjAP5NSunXzeyPAXwSwBMAHgTwZOe+5jXqN3+yu6DwgpWxyEXiWTH0iTKuORILilGw0BKioq2t5ccd+bk0R9fzdUiAa0ZRMBAdh6o9l29794orNbP4V68Fcxux0EeBPz0Simbr/mVyvkYJRHzoKGaqI2AoqkDE27CgyIlLQHfl5qhyllu2RE5bVfgd8zJOTIoCit5XvZ2Nbyk23TrMifJy5zoR1xJn8FkA/9rMXsWOhvCla9iXEOKA2VfcYkrpOQDPtf//AwAfuv5TEkIcBIpAFEIAOIAuzIuwPgB4f79Yz31wlwiEIPiH14kScDj4J0qA4gq+6/lxmiP+68h8Y0Dj3E9vBoFfzt2RqCtTGNpJmkE9pmShoT9nTv6Zj0nPCKQV18Vo5NeZbexdtThKIHLBPz26a7Fm0HA8Vw/NwHVLijoqcQVlvi5RcRPSA3gMABvUZXm9zMenqnfcNqfKSzTOt4lKzRwPApH6oDcDIQQAGQMhRIuMgRACgIyBEKJFxkAIAUDGQAjRImMghAAgYyCEaFltdeSqRHnLid3xkQ23ThpRiA0FJjVRAs6AOwVx0JGfS7FFrdKj7k6zPFqGk4zmR3w40OxYPpfZGnUK6lGtp5pQRaUop4qCjrgiUVRd2M2FYqaiufVpYc77CZOMCKPL3SfoiAOI+LjNsEcAUY+5JQpE4rFr0Q6g4m5Jpa9OdbTKE85ura5k49uCBCMOMjpZ5M/ceuGfwdtTnwguj94MhBAAZAyEEC0yBkIIAKtOVCpL2PFj7w3r2466VZph7nPP1yjRZ+TtFyf/cOedqMBFtZ1rD8Us8ANJV5iPOenIz2W6Qck/g+5EH5dMQ8lBxazbZrvjBAlRM5JouIqxK0qCfkk6fE4p8Kndfrm4SdOj8jEJDfN11k2CbXguSxRYToPcBy/GXg8YD/Jlx4JOyNwxKaps3EVDlWJmyQtKxTInCb0ZCCFaZAyEEABkDIQQLTIGQggAqxYQiwLN0d0Kw5MTXk2rKTDGVeMJZsyC4bxHewbuZFVte9FrMN47+GceBPa4ajyusk73NjULZbNuQagpu8VAJ7TSOqG46dqR++tUc9wLd4uLgoHcKe2/dX1aoxsSzM0qalVH1zbV3X8PyyEFFI18+7Ijg+5W6cwb81xAj1qyM7cUW9n4BI0BYCOoct0HvRkIIQDIGAghWmQMhBAAVqwZpMJQb+w6mPN1X453ts7+PwUQBTPmdTiBJaxmS65WE5QgtmZvWxkm7ZBfzj53WMGX1ilII+BgoZ250X7Zty/9SfM5smYQXYNmwEk7fp16ROsMaXKBLx8l+3RC/n45yn3ssvQJOoPB3n54HWgGBVU2Hlb5PsZDrxmsVfkiWJaUAAAFMUlEQVSytdKvM6MbvUkizY+ak26bK1QCmgOVouSmPl2XIvRmIIQAIGMghGiRMRBCAFh1nIHlXYnDmIGub+HBt33WCDiBJUpU8t/l/Uol5Zqwnx5pBl4j6NY8eJu58+WDE+D8G+5iHHVHYl2BC39E3YVom3rD++WcyANO7OHfAyiK/Fj8/d+Kbk1hTN/72bcHgFHlk4q64G7JXKiE9QEAODXu7nxc0DmenxzPj1P4uf7N9Fg2Zi3iOItfAE4OLtGS/9M5N0BvBkKIFhkDIQQAGQMhRIuMgRACwKqDjsxQLwiIYWBMR/BMFBjDiTIucSbM2yDBKlAZ3bF6FJ31iUk0jhKIWOxziT5+GxYM+7Q45+O4IKQooGhMYmwULNQhGFZB4I8TDHvk1nBQEQuGkVjIy1gMbIL7vk6t04dU+TgKKBoU+TpVkHR0ZT6icX4DXp/7G110lI0+HlRUujDwFcT6oDcDIQQAGQMhRIuMgRACAGApLZEwsuzBzN4A8CMAJwG8ubIDXxs301yBm2u+N9NcgZtrvotz/TsppVNdG6zUGLx3ULOzKaUzKz/wEtxMcwVurvneTHMFbq75LjNXuQlCCAAyBkKIloMyBo8d0HGX4WaaK3Bzzfdmmitwc81333M9EM1ACHH4kJsghACwYmNgZg+Y2ffM7FUze2SVx+6DmX3ZzC6Y2YsLy06Y2TNm9kr789aDnOO7mNndZvZVM3vZzF4ys0+3yw/rfMdm9jUz+3Y7399pl7/fzJ5v5/sVMwuCrw8GMyvN7Ftm9nQ7Psxz/aGZ/YWZvWBmZ9tl+3oWVmYMzKwE8J8A/FMAPwfgU2b2c6s6fk9+H8ADtOwRAM+mlO4F8Gw7PgzMAfxWSukDAD4M4F+21/OwzncC4P6U0s8DuA/AA2b2YQCfB/CFdr4/BfDQAc6R+TSAlxfGh3muAPCPU0r3LXxS3N+zkFJayX8A/iGAP18YPwrg0VUdfx/zvAfAiwvj7wG4s/3/OwF876DneJV5PwngYzfDfAGsA/gmgH+AncCYKnpGDniOp9t/QPcDeBo76W6Hcq7tfH4I4CQt29ezsEo34S4AP14Yn2uXHXbuSCm9BgDtz9sPeD4OM7sHwAcBPI9DPN/2tfsFABcAPAPg+wAuppTeTSU8TM/EFwH8NnZzVW/D4Z0rsJOG+9/N7Btm9nC7bF/PwipTmKMkVX3KuEbM7AiAPwHwmZTSO9YnF/iASCnVAO4zs1sA/DcAH4hWW+2sPGb2SwAupJS+YWYffXdxsOqBz3WBj6SUzpvZ7QCeMbPv7ncHq3wzOAfg7oXxaQDnV3j8ZXndzO4EgPbnhQOez3uY2QA7huAPUkp/2i4+tPN9l5TSRQDPYUfruMXM3v2jdFieiY8A+GUz+yGAJ7DjKnwRh3OuAICU0vn25wXsGNoPYZ/PwiqNwdcB3NsqskMAvwrgqRUef1meAvBg+/8PYsc3P3Bs5xXgSwBeTin97sKvDut8T7VvBDCzNQC/iB1x7qsAPtmudijmm1J6NKV0OqV0D3ae0/+RUvp1HMK5AoCZbZjZ0Xf/H8A/AfAi9vssrFjk+DiAv8KOr/jvDlp0Ceb3hwBeAzDDzpvMQ9jxFZ8F8Er788RBz7Od6z/CzmvqdwC80P738UM8378H4FvtfF8E8O/b5T8D4GsAXgXwxwBGBz1XmvdHATx9mOfazuvb7X8vvftva7/PgiIQhRAAFIEohGiRMRBCAJAxEEK0yBgIIQDIGAghWmQMhBAAZAyEEC0yBkIIAMD/BT9zl4h7oN6OAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data_raster = banddata_raster.ReadAsArray(region_x_offset, region_y_offset, region_band.XSize, region_band.YSize)\n",
    "\n",
    "print(f\"{data_raster}\")\n",
    "plt.imshow(data_raster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-- -- -- ... -- -- --]\n",
      " [-- -- -- ... -- -- --]\n",
      " [-- -- -- ... -- -- --]\n",
      " ...\n",
      " [-- -- -- ... -- -- --]\n",
      " [-- -- -- ... -- -- --]\n",
      " [-- -- -- ... -- -- --]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fcfbb4879b0>"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQMAAAD8CAYAAABzYsGzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHa9JREFUeJztnWuMZdV153/rPqqqq7uhbSiTDo3cJEGWIzkBqcfB4/lASIgIE9lR4iQ41oiJUFCURLEdWzYkM57xJErsKDFWpDwGxZb5YBk7L4FQogQRk4cUYTcGBwjGTQiOWzB026aBpqvqvvZ8qNO49tqLOufeqltVt/n/pFb13rXP2euee2vds/9nrbUtpYQQQrR22gAhxO5AzkAIAcgZCCEq5AyEEICcgRCiQs5ACAHIGQghKuQMhBDAJp2BmV1rZo+b2RNmdvNWGSWE2H5s0ghEM2sDXwWuAY4DXwTemVL611c65sILL0yHDx+eaL5zhWOnv170dWzk2sOsfcnid03VJnFu88ADD3wjpbRUN66ziTneDDyRUnoSwMzuAN4OvKIzOHz4MEePHt3ElLPPdf/wK0XfgbnlrH3R/AtZ+9YrPjtVm8S5jZl9rcm4zSwTLgbWf80dr/qEEDPIZpyBBX3FmsPMbjKzo2Z29OTJk5uYTggxTTazTDgOXLKufQh42g9KKd0G3AZw5MiRV32K5HMre2rHzLcGWftjj11TjNnfWsnaXcuP6afyrf35N/xjExPFq5TN3Bl8EbjMzC41szngeuCurTFLCLHdTHxnkFIamNkvA38DtIFPppQe3TLLhBDbymaWCaSU/gr4qy2yRQixgygCUQgBbPLOQIxPf9gu+kYpfzCzOsrflq4LQlrrywXDtnltdoDnU199S25LIDJ69rfzGIjrv+fVHSdyLqM7AyEEIGcghKiQMxBCANIMps6b7vpQ1t63UK7/O608Uem87koxpo6XRvNZe3XUHfsc57dfKvraLqj02NcPFmMuaucfo32thazdT+Vrnj/45Nj2iemiOwMhBCBnIISokDMQQgDSDKbOvoXVrP2aheVizAXz+Vp9ae7FrB3FGbw4yhOe+imPXzg9zNftW0Uvld8fK04TaI16WXs1lTEPzxzPtYdvDXON4/8NzyuO+Y/+BVn7F97w9xsbK8ZCdwZCCEDOQAhRIWcghADkDIQQFRIQp8wl+09l7fO7pYD4nfPPZ+2FVj9rRwFEQ1d17rnB3qx9epAHIQGcHpZ96znR2l/07WvnAugFndPFmNWU2z/vqy4F3zlnXJDUqdFi1v7mcF9wzFzW/rkv/FzWfnaltP+by/l5T6/k8z769g8Xx7xa0Z2BEAKQMxBCVMgZCCEAaQZbjl/H7nW1TA4EmkHb7ajUcm2/VgZ4fpAHHZ0a5GvjF/pl0NGZQa49DFyg0lyrDA7a3801gyfnLirGvNTN1+F7W/kxURGVoQteOjXM7V9JpU7iA6mWXaDS6rCcZ7Wf9/V6efvST/9WccyoXxagWU+rWwaB/fu7fm3DY2YB3RkIIQA5AyFEhZyBEAKQMxBCVEhAFGIL8BWt9rtsVYDv2Jvvrn1oMQ9I+/0rPrP1ho2B7gyEEICcgRCiQs5ACAFIM9hyDnTPbPj7waj0v8NW3vdc3yUdBQlGLw7yAJwX+/mYM4MyUGl5UFcxuZxnxQX2HOu8rhjzLZck5ROt5l0byspMvppTmJzlApVabhepFn5XqXpGg+D7cJAngdHKz5tG7vdA21W43tvtFWN8AFeb/JjffvS64hh/7fx1WurkOgRMvuuV7gyEEICcgRCiQs5ACAFIM9gUN3/5J4u++Qnc63MuyajndmH+Vi//PcBppxG81M81gl6w23NvsPHbPQzWwivumIX2+cWYl1whlf1uR6j5IAHK6wp+vT8isGXkE63yix0dk9wO18Ohe4OG5TH46+CHBIe0na7gNQSA/ih/T4buu9gnnwHMt/LX7K9bVDn7vqcuKw1sgO4MhBCAnIEQoqLWGZjZJ83shJk9sq7vtWZ2j5kdq36+ZrpmCiGmTZM7g08B17q+m4F7U0qXAfdWbSHEDFMrIKaU/sHMDrvutwNXVf+/HbgP+OAW2jUTRMJYVJVoPaf6pUj0oqtK1HNC0zdX8qAeKIU9X9GnHwiInkJMi8a4IKnn50r7PV7IiyooeWHMC2FeXAMYpUC52yyBaEqh/eVjUiA6evHVi4VQCp7Pu89Cr13+OXbcddnTrg/ouurwMdfT7LpNqhlclFJ6BqD6WYalCSFmiqkLiGZ2k5kdNbOjJ0+enPZ0QogJmdQZPGtmBwGqnydeaWBK6baU0pGU0pGlpaUJpxNCTJtJg47uAm4APlL9vHPLLJohoiAXn4j0LZd09HyvrFrsk4HOuACi55fLY7xGMIgSbhyjIElqXJb7ZQJRt5Wva/3afq5dBsZ0XAVor790WuUxA7cOj5K+tgJz9iefABVqBrkt0XV6oZW/j4OOO2ZYHuODsQ7M5dW1F1tlQNqkNHm0+Bngn4E3mNlxM7uRNSdwjZkdA66p2kKIGabJ04R3vsKvfmiLbRFC7CCKQBRCAEpUGov3PvgzWXtIfcyALzry3Gq5xvMxA2dWc81gebmMXRg6jcAn5ETPwutoz5fr9CLRJ3jW79fLPXN6RiqfufvYA/8MntKUiYqXFOdwBVGCaQps6DWE8ju05zScM0GikscXm+kExxR6TIMYjknRnYEQApAzEEJUyBkIIQA5AyFEhQTEMXh29bys7cUdKAXDF12Q0Uu9Ugz0gmFv1SUdLQdVjb1AOIFgiBcMA3FwFCXyOPxW6F7kioS/gasI7YW9CD/GBzc1CkLy8zS5bE7X84IiQG8lvwYWvB4vtPpqSJ12ICC6gC1/js+95f+W9k6I7gyEEICcgRCiQs5ACAFIMxiLU736wh6+qMhpl3Tk9QGAlZVcExiuuiCdldJnm9/1pwGps3HCTaqPk2HQoGhK361rI0t9gE20xq47xmsRPkAKoNvJ19ydTn6OfrucN7m+Vs9d/8DUkdMMyj2YYdDJr535isqBZuD7mhStmRTdGQghADkDIUSFnIEQApBmsCGX3/0/svZig2Kgfjm53Mv1gF6vXPON+s4nu3b7TOCz3frYL7lHXh+AYvHudxOO1tx+TG8QFPp0hVX9+t/vNgTQK8Y4PSDQEPx62T+Db4Jfg1uwTsfN4zctsiD2Ii273Z1S+ac16rq53XUZBNfJnH7RpIjtpOjOQAgByBkIISrkDIQQgJyBEKJCAqIQM8RXfuJ/Te3cujMQQgByBkKICjkDIQQgzWBDeq5qsQ/KaZJc488R7WqUfNDOch700l6t3ynYFyBuBelBRY6LDzoKinaMOm534SDoKEqwyY6JKh27a+fb0bX1xT/8iGiX5ibBTHX4pDALXm7bJZNFVaSTO48vCJ2CoCOivimhOwMhBCBnIISokDMQQgDSDF7m9beXe8d2F/LLM+yM7zuHbh1eJCUBuOIZ7WVXUDTSDHzSUXIJLcE0fq1bFDOJCqK6nZt8IVMok2daE6xzvUYQnWM4yg32CVJRQlR5jgbvoVvb+7q3rV5wnYZ+7mhMq25IQVGQZorozkAIAcgZCCEq5AyEEICcgRCiQgJiha8oAzDotzdsh+dxQlhRTSiodGROkGr1899Hu257rS/Y9Tw4yM3rt3H3VYApRcahBfbX7aY+wW5J0THm5q4LdgIYOAF06N7DtFz+CXjBsAg6Ct6P9tBXngrs72xcnSrQb6E1wU5ZE6I7AyEEIGcghKiodQZmdomZfd7MHjOzR83s3VX/a83sHjM7Vv18zfTNFUJMiyaawQB4X0rpS2a2H3jAzO4B/jtwb0rpI2Z2M3Az8MHpmTpdoiQdBjW+MlgKh+u+9fTLAT6IpeXGRJrByC3d/bq2FWggPpimKPIbfBqS+75oErPTRCPwlYG9LVGiUqFNuKSpUbDLlE8C8wFeraDydOel/Dxttz2S13TWjPMTB++zO08TzeArv/HeYLLpUPvWppSeSSl9qfr/i8BjwMXA24Hbq2G3Az8+LSOFENNnLM3AzA4DVwD3AxellJ6BNYcBvG6rjRNCbB+NnYGZ7QP+HHhPSumFMY67ycyOmtnRkydPTmKjEGIbaOQMzKzLmiP4dErpL6ruZ83sYPX7g8CJ6NiU0m0ppSMppSNLS0tbYbMQYgrUCohmZsAngMdSSh9b96u7gBuAj1Q/75yKhVPi9Z/6aN4RCX8+kMSLf5FOVqsglrS8+OcEqijIpUjQ89NaZMfGB43C1+ObEzyNDrIJ6y5TalKRyAu8QQyS1byHneXSEF9ZqtXLfx8KiIw/JrJ3J2nyNOGtwH8DHjazh6q+X2PNCXzOzG4E/gP4qemYKITYDmqdQUrpn3jlzOsf2lpzhBA7hSIQhRCAEpW+TbCItVUXoOKCg8Jl7QTrwGJNOti4DZDG342cIDImb0ZbjbvXE1Ugrg0yivSLCb6GvI7QchWJvfYCFO+Hr/bkq0pBGWTUWXHnGAQBUf5StuurUxV6zA7/NerOQAgByBkIISrkDIQQgDSDbxOs9YsiI/75/0Tr9pIiyWjgfx8c4zu8Ww+WrH65b65wRisInBi5E0UFiH2RlGLeKM6ggb1183RPOw0nSjbzNVPce9Z2eg1AZ9mNWXG7PQXve1FcpkGxmVE3bz986/YlJUXozkAIAcgZCCEq5AyEEICcgRCi4lUhIL7+T36n7PRuMAo6qhX26lWvYkuzQFgqgop8YEw/KqlUzFRrSx1eLIRSVIzCi+qE1BQE4EyyaZivWtxe2TjBC0rbfHBQdEx7OR9UiIxRUaw51xF8zQ4W83Yhou4wu8wcIcROIWcghADkDIQQFXIGQghAzkAIUSFnIIQA5AyEEBWvijiDqGhH8aw4SHLxMQJFcZPg+XrtLjlRoo+fx523UXHNYtff4Nm+62r5r4JgHr+L9KgTFPZwY4rYiugrp64gSrQjkbOv7YqO+KIkAK3BxklG0XvY7o0fBTHq5i+yv68cM9jjjpmfJNpieujOQAgByBkIISrkDIQQgJyBEKLinBQQD//B7+Yd3WCQ126iSkc+MckLfU121vFiWnTFfQVfJ3q1okQlhxfxUiD00Rs/mcnvstRuILQWFZXCaTfencqfE8r3o3PGX6fymHZRncoJihNUs+4vli+od17eHuwtr/9gT9737+953/iTTxHdGQghADkDIUSFnIEQAjhHNYMmRUd8kZHwmLRxcZNojRrtfrSeURDkUhdkVHdOgKFz6+HOzTW7MEdxQOY/IeNv7hxTU7U4Cgby18nvdBRpK8V7NqzXX5KrGj1YyNvD+fIF9vfl5x2cX4oRac8WldOeErozEEIAcgZCiAo5AyEEcK5qBsX6M0jaaTdIOvLnKdaf0THu2bePVWiwFvbPxn2yTTPK1+ynNr+rcXCMtzfcPMl/pdQkXq2dyM/jbImOceftLPtrHWgGNRqB1wcAhq64qU8wGuwtz+PjCiJ9oL3QQPzZQXRnIIQA5AyEEBW1zsDMFszsC2b2ZTN71Mw+XPVfamb3m9kxM/usmfnK8UKIGaLJncEqcHVK6fuBy4FrzexK4KPArSmly4DngBunZ6YQYtrUCogppQScrprd6l8CrgZ+tuq/HfjfwB9tvYn1XPbbH8vanWhrbkcRdNRAQGyyA04R5OI1o0DPKoKOnBDWXm2QqNREY/QZQ645iozzQmtYqalG/IuO8a+5EH2jikp5u7Oad0RBR8m9Zv++D7vlZ2Uwv3GQkRcYAdJ8bou1y6CjJ6//9fLAXUQjzcDM2mb2EHACuAf4N+BUSunsR/04cPF0TBRCbAeNnEFKaZhSuhw4BLwZeGM0LDrWzG4ys6NmdvTkyZOTWyqEmCpjPU1IKZ0C7gOuBA6YvRy5fgh4+hWOuS2ldCSldGRpaWkztgohpkitZmBmS0A/pXTKzPYAP8yaePh54B3AHcANwJ3TNHQj2ssbV+eNClj49X+kB9QVM/FFSKA+yChYShZrXV+dNww6Kqoh1/v1cgfoBhlFNdWFoV4jiAJ/Ss2gvuiI1ybaK3l71ImCy/L2YCG/ToM9gWawN+/ruyCjJlWNW1FxmV1OkwjEg8DtZtZm7RP3uZTS3Wb2r8AdZvabwIPAJ6ZopxBiyjR5mvAvwBVB/5Os6QdCiHMARSAKIQA5AyFExTmRtdju5e2W22YrCpQpsu+aBBTVBMpAKQZ6IawdVuPZWGyyQammla/JjylfUNvHHBXBTsHkTaoW14h/XviLzltcgwZVi30wUCQg+ipFXjAcBJWO/dZoI1dde7gQGNOaPcHQozsDIQQgZyCEqJAzEEIAM6gZvOm9txZ9Hb9V9wRbakfxNyMXsOLX1F6rAOi4pCKvIbQC2+q2cW/3gwW0X6e78shhoo8bM5obf4eleBv6jdf7YXKT00Fqt7IPxgwW849vb3/53TZw6/vhno0DigBGcy6YySUmRdvSe82g1d7dlZAjdGcghADkDIQQFXIGQghgBjUDIWaBJ376f+60CWOjOwMhBCBnIISokDMQQgByBkKIipkTELuny4CP2u3Iwoq+bkjkFn3giztPFNzkKxm3V/LgkyYBRD6jqLVcZhD55J/UyV9Aq12+oJG7LqlXH3Tk5wmTjhy+InERlAS1rzlisNB2bb91enlMf79PVMp/P+qWtvnEJF/ZaDRfvoftXb7dehN0ZyCEAOQMhBAVcgZCCGAGNIMr3/l7WXshSMDxGkGxs05UabdBLpMvluHXvlHSTns172w5zaDVC9b/7jX57eJby04UgWJ3JL9Op1P6+Zbri7YjLw9y1yAotFJ7nuD3hb3+LZsr7ffzeM3AFzuBcvej5JKMktMH1uZ2Y1ybQGfozu3u7daboDsDIQQgZyCEqJAzEEIAM6AZLHwzXy+Hu/P0Ny6UQfBs3IYNKm62Nl6XR+fwfa0Vt5ZcLSuiFDEDfp3eD9ajLVfMxGsInWL7pGIMPhbB/z6gWOsDzJdzZccE3znmukZOzxgFcQcjrwn4XbEiM7w04dvB16E/T2r7mIjy8/T4T34omHy20J2BEAKQMxBCVMgZCCEAOQMhRMWuFxDnv7Gcd4S7C/myuS4aKEqUGTUREGsEtuC8Vjf3itvuKcDcvGlQH9BinfyttH7g552omIJkpgJ/DeaCKB1H6tbPk1xA17CbjxkulGrg0AUiDbtONN1Yx9w6mkSszSC6MxBCAHIGQogKOQMhBLALNYMf/e73Z+2WX4MH6+fk1/8NCnCQGmgG7XwR6tfyIW7t3kibcBQawbC+cIZ/xV5DWLPFBTP5QKUo4cjrDAvleYd7cx3BBxBFuyONnEbgk7OGC0HQkZ+6ZldpKHfBKl5jsHuyD0Qa+eCzJlrLDHJuviohxNjIGQghgDGcgZm1zexBM7u7al9qZveb2TEz+6yZzdWdQwixexlHM3g38BhwXtX+KHBrSukOM/tj4EbgjzZrUHruVN72STpRzEC75gHzBOt2AOvmlye57BpbmK89R9qTj4kKihaJSX5MlKjkdQX//D+4Jsm/noWu+315zMj1DfYGmoGLCSiKjASagX/Nfsfr6JjhBLtGe7yuEO4q7bQVXxznazd+YNN27EYa3RmY2SHgvwJ/UrUNuBr4s2rI7cCPT8NAIcT20HSZ8HHgA3zbR14AnEopnf16Og5cHB1oZjeZ2VEzO3ry5MlNGSuEmB61zsDMfgw4kVJ6YH13MDR8npdSui2ldCSldGRpaWlCM4UQ06aJZvBW4G1mdh2wwJpm8HHggJl1qruDQ8DT0zNTCDFtap1BSukW4BYAM7sKeH9K6V1m9qfAO4A7gBuAO7fCoOELp/OOUa7wWLd8aFHcpnjxLAoWmkBUtEW3Hc98acvovMW8PdfesA0UuwmZq/bcOeWStaAQUr34N9pTJhQNXcDQcL6+WvLIiXaDhfJa9hedGOguSyTSFW+afzuie0+v6xWVj8tDooCn3LZyQLmL1MbnOFfYTJzBB4FfNbMnWNMQPrE1JgkhdoKxwpFTSvcB91X/fxJ489abJITYCRSBKIQAdjhR6ZrO9UWfD/SxuXwNHiXgFME/fkyUgOMDe6ICIu48aTHf5ne0v9z2d7AvX8gO9rq1fLe0xa/V26sNFqnO/lER+BMUB3HJP/09XjMopykSiII4q/6+jQOIWpFmMGowxtvi7PO2FIlMwTF+R6VGxwTVkM9FdGcghADkDIQQFXIGQghAzkAIUSFnIIQA5AyEEBVyBkIIQM5ACFGxo0FHnYsPFn1p3iXYuAo+UQJOkaTjqvP6RCCA1nIeZBRur76aR8KMFvO5fVVggN75+SXtL7pdgBpU62mv+q3fg0EuUclXG46qC3tbBi7vKgqIKioQBV8fAxd7VQQHlYcUmxI1SQbySUfDPe4aBIX3igCioBqyZ+Te1qd++f3xwHMM3RkIIQA5AyFEhZyBEALYZs3g2MPHsx2ThhcdKMaM5nOTBnucHjBf+i+/1h25HX6joh3t1a5rl2vJjtMVBgs+6ai0pbd3Y40gSvSpqxTc7tX7bD9PtP7v+YQin+gTbLDcqGrxgiu00mA3ZL/+D3URj7sMg0U371xQefpVkmS0FejOQAgByBkIISrkDIQQgJyBEKJiWwXEy950iL8++rubOsd/uuFjRZ+vVjPYM35gT3e5FJq6Z9x53BBfXRhKMdBX7E2dqNKRO6/bKt0HC0V40S4UA7sbt33wUHReXykIYOirFLd95eZAxNuCr6G0x1XO7pSRS23X13LRTqORvg/PoishhADkDIQQFXIGQghghxOVJuGLt//qTpswFpf/Yq5xNAnsafXzdn9veYxP9Bn5ZKxoEym3th8WGkK5ti90hiiwZ95rBG7t7tuATRAMZO5Fd+dyzaDbLSOX5jpB1et1DAPN4OG3/Z+xbTsX0J2BEAKQMxBCVMgZCCGAGdQMZo2H/nB2NI7v+Z0yhqPYgWgxqEIy557lu7V7O1jLt51m4PWAVoMiJIvzvdyMdjnP3m6v6BMxujMQQgByBkKICjkDIQQgZyCEqJCAKF7miQ/MjtgJ8J//9oNZOxIL513Q0WInHzNQotLL6EoIIQA5AyFEhZyBEAIAS2n7qsea2Unga8CFwDe2beLNMUu2wmzZO0u2wmzZu97W16eUluoO2FZn8PKkZkdTSke2feIJmCVbYbbsnSVbYbbsncRWLROEEICcgRCiYqecwW07NO8kzJKtMFv2zpKtMFv2jm3rjmgGQojdh5YJQghgm52BmV1rZo+b2RNmdvN2zt0EM/ukmZ0ws0fW9b3WzO4xs2PVz9fspI1nMbNLzOzzZvaYmT1qZu+u+nervQtm9gUz+3Jl74er/kvN7P7K3s+a2VzdubYLM2ub2YNmdnfV3s22PmVmD5vZQ2Z2tOob67Owbc7AzNrAHwA/Cnwv8E4z+97tmr8hnwKudX03A/emlC4D7q3au4EB8L6U0huBK4Ffqq7nbrV3Fbg6pfT9wOXAtWZ2JfBR4NbK3ueAG3fQRs+7gcfWtXezrQA/mFK6fN0jxfE+CymlbfkHvAX4m3XtW4Bbtmv+Mew8DDyyrv04cLD6/0Hg8Z228RXsvhO4ZhbsBRaBLwE/wFpgTCf6jOywjYeqP6CrgbtZ24x+V9pa2fMUcKHrG+uzsJ3LhIuBr69rH6/6djsXpZSeAah+vm6H7Skws8PAFcD97GJ7q9vuh4ATwD3AvwGnUkpnUwt302fi48AHgLM13S5g99oKa5v//a2ZPWBmN1V9Y30WtjOFOdoAUY8yNomZ7QP+HHhPSukFs/p9JneKlNIQuNzMDgB/CbwxGra9VpWY2Y8BJ1JKD5jZVWe7g6E7bus63ppSetrMXgfcY2ZfGfcE23lncBy4ZF37EPD0Ns4/Kc+a2UGA6ueJHbbnZcysy5oj+HRK6S+q7l1r71lSSqeA+1jTOg6Y2dkvpd3ymXgr8DYzewq4g7WlwsfZnbYCkFJ6uvp5gjVH+2bG/CxspzP4InBZpcjOAdcDd23j/JNyF3BD9f8bWFub7zi2dgvwCeCxlNL6ssa71d6l6o4AM9sD/DBr4tzngXdUw3aFvSmlW1JKh1JKh1n7nP5dSuld7EJbAcxsr5ntP/t/4EeARxj3s7DNIsd1wFdZWyv++k6LLoF9nwGeAfqs3cncyNpa8V7gWPXztTttZ2Xrf2HtNvVfgIeqf9ftYnu/D3iwsvcR4ENV/3cBXwCeAP4UmN9pW53dVwF372ZbK7u+XP179Ozf1rifBUUgCiEARSAKISrkDIQQgJyBEKJCzkAIAcgZCCEq5AyEEICcgRCiQs5ACAHA/wd1yB8aYsdnWAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Mask zone of raster\n",
    "zone_raster = np.ma.masked_array(data_raster, [x == region_band.GetNoDataValue() for x in region_array])\n",
    "\n",
    "print(f\"{zone_raster}\")\n",
    "plt.imshow(zone_raster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1187, 1368)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.min(zone_raster), np.max(zone_raster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 81, 144, 258, 386, 403, 389, 371, 309, 144,  65]),\n",
       " array([1185., 1204., 1223., 1242., 1261., 1280., 1299., 1318., 1337.,\n",
       "        1356., 1375.]))"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.histogram(zone_raster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1269.4353701527614,\n",
       " 1269.4353701527614,\n",
       " 1269.0,\n",
       " 38.35120031548615,\n",
       " 1470.8145656385454)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.average(zone_raster), np.mean(zone_raster), np.ma.median(zone_raster), np.std(zone_raster), np.var(zone_raster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# raster to vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdal.UseExceptions()\n",
    "\n",
    "src_ds = gdal.Open(\"region.tif\")\n",
    "srcband = src_ds.GetRasterBand(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "srs = osr.SpatialReference()\n",
    "srs.ImportFromWkt(src_ds.GetProjection())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "dst_layername = \"POLYGONIZED_STUFF\"\n",
    "drv = ogr.GetDriverByName( \"ESRI Shapefile\" )\n",
    "\n",
    "if os.path.exists(dst_layername+\".shp\"):\n",
    "    drv.DeleteDataSource(dst_layername+\".shp\")\n",
    "\n",
    "dst_ds = drv.CreateDataSource( dst_layername + \".shp\" )\n",
    "dst_layer = dst_ds.CreateLayer(dst_layername, geom_type=ogr.wkbPolygon, srs = srs )\n",
    "\n",
    "df = ogr.FieldDefn( 'id', ogr.OFTInteger )\n",
    "dst_layer.CreateField( df )\n",
    "\n",
    "maskband = srcband.GetMaskBand()\n",
    "# print(maskband.ReadAsArray(0, 0, maskband.XSize, maskband.YSize))\n",
    "\n",
    "result = gdal.Polygonize( srcband, maskband, dst_layer, 0, [], callback=None )\n",
    "\n",
    "dst_ds.Destroy() # 注意关闭数据源来写入数据\n",
    "src_ds = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shp = ogr.Open(\"OUTPUT.shp\")\n",
    "lyr = shp.GetLayer()\n",
    "featList = range(lyr.GetFeatureCount())\n",
    "\n",
    "for FID in featList:\n",
    "    feat = lyr.GetFeature(FID)\n",
    "    geom = feat.GetGeometryRef()\n",
    "    print(geom.ExportToWkt())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# use rasterio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "from rasterio.features import shapes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = rasterio.open('region.tif')\n",
    "dataset.indexes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "band1 = dataset.read(1)\n",
    "band1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shapes(band1, mask=None, transform=dataset.transform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.nodata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[x == dataset.nodata for x in band1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open('region.tif') as src:\n",
    "    image = src.read(1)\n",
    "    results = []\n",
    "\n",
    "    for i, (s, v) in enumerate(shapes(image, mask=None, transform=src.transform)):\n",
    "        if v != dataset.nodata:\n",
    "            results.append({'properties': {'raster_val': v}, 'geometry': s})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = [x['properties']['raster_val'] != dataset.nodata for x in region ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lst = []\n",
    "for x in region:\n",
    "    if x['properties']['raster_val'] != dataset.nodata:\n",
    "        lst.append(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(lst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "src_filename = 'region.tif'\n",
    "dst_filename = 'out.shp'\n",
    "\n",
    "result = subprocess.run([\n",
    "    'gdal_polygonize.py', src_filename, dst_filename, '-b', '1', '-f', 'ESRI Shapefile', '-q', 'band', 'id'\n",
    "    ], stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "result.stderr.decode('utf-8') == ''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
